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Abstract 

This paper is a contribution to the old problem of representing a signal in the coordinates of 
time and frequency. As the starting point, we abandon Gabor’s complex extension and re¬ 
evaluate fundamental principles of time-frequency analysis. We provide a multicomponent 
model of a signal that enables rigorous dehnition of instantaneous frequency on a per- 
component basis. Within our framework, we have shifted all uncertainty of the latent 
signal to its quadrature. In this approach, uncertainty is not a fundamental limitation of 
analysis, but rather a manifestation of the limited view of the observer. With the appropriate 
assumptions made on the signal model, the instantaneous amplitude and instantaneous 
frequency can be obtained exactly, hence exact representation of a signal in the coordinates 
of time and frequency can be achieved. However, uncertainty now arises in obtaining the 
correct assumptions, i.e. how to correctly choose the quadrature of the components. 

Keywords: Hilbert Space, Signal Analysis, Instantaneous Frequency, Hilbert Spectrum, 
Latent Signal Analysis, AM-FM Modeling 


Part I. Latent Signal Analysis and the Analytic Signal 

Interest in the proper dehnition of instantaneous frequency hrst arose with the 
advent of frequency modulation for radio transmission in the 1920’s .... [T]he 
“analytic signal procedure” devised by Gabor results in a complex signal that 
has spectrum identical to that of the real signal for positive frequencies and zero 
for the negative frequencies .... However, instantaneous frequency is a primitive 
concept and not a question of mere mathematical dehnition .... One should keep 
an open mind regarding the proper dehnition of the complex signal, that is, 
the appropriate way to dehne phase, amplitude, and instantaneous frequency. 
Probably the last word on the subject has not yet been said.—Gohen [1] 
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In the first part of this paper, we present the Latent Signal Analysis (LSA) problem 
as a recasting of the classic complex extension problem. Almost universally, the solution 
approach has been to use the Hilbert Transform (HT) to construct Gabor’s Analytic Signal 
(AS). This approach relies on Harmonic Correspondence (HC), which may lead to incorrect 
Instantaneous Amplitude (lA) and Instantaneous Frequency (IF) parameters. We show that 
by relaxing HC, the resulting complex extension can still be an analytic function and we can 
arrive at alternate lA/IF parameterizations which may be more accurate at describing the 
latent signal. Although the existence of other lA/IF parameterizations is not new, Vakman 
argued that the AS is the only physically-justihable complex extension. We argue that by 
modifying the differential equation for simple harmonic motion [2, 3], our parameterizations 
are also physically justihed. 


1. Introduction 


Many physical phenomenon are characterized by a complex signal 

z{t) = x{t)+jy{t) = p(f)e^®(‘^ 

d 


( 1 ) 


where p{t) is the signal’s lA, Q(t) the signal’s phase, and f2(t) = —0(t) the signal’s IF. We 

dt 

will call z{t) the latent signal because only the real part, x{t) is observed and the imaginary 
part, y(t) is hidden, i.e. the act of observation corresponds to the real operator 


x{t) = iR{ 2 :(f)}. 


( 2 ) 


It is often desirable to analyze the latent signal—which is known to completely param¬ 
eterize the physical phenomenon. Thus the complex extension problem becomes that of 
determining z{t) from the observation x{t), i.e. recovering the quadrature y{t). In the clas¬ 
sical approach, we seek a unique rule /!{■} such that the estimate of yit), y{t) = C{x{t)}. 
These relations and the complex extension problem are illustrated in Fig. 1(a), where we see 
many latent signals mapping to a single observation under the real operator. In the LSA 
problem, given x{t) we seek z{t). 

In the context of time-frequency signal analysis, we desire the instantaneous parameters, 
p{t) and Q{t). Once we have determined a rule for y{t), the instantaneous estimates are 
given by 

p{t) = ±\z{t)\ = ±^/x^{)~+WW ( 3 ) 


and 


n{t) 


d 

dt 


arctan 



( 4 ) 


The very dehnition of lA and IF hinges on y{t) and hence /!{■}. 

The complex extension problem is well-known. In 1937, Carson and Fry formally dehned 
the IF based on the phase derivative of a complex FM signal [4, 5]. In Gabor’s seminal 1946 
paper [6], he introduced a Quadrature Method (QM) as a practical approach for obtaining 
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(a) 


(b) 


Figure 1: Set diagrams for the LSA problem, (a) There are an infinite number of latent signals, z{t) that map 
to the observed signal x{t), according to x{t) = 3fi{2;(t)}. We seek a rule, £{•} to determine an appropriate 
latent signal z{t). (b) Almost universally the rule chosen is thus this limits us to only a subset of 

latent signals as illustrated by the shaded part. 


the complex extension of a real signal and showed its eqnivalence to the HT. In this context, 
the rnle is given by 

y{t) = (5) 

where 'H{-} is the HT operator and 

z{t) = x{t) + jV.{x{t)] (6) 


is termed the AS. This is illnstrated in Fig. 1(b) where we see nsing the HT to estimate the 
qnadratnre leads ns to a snbset of the latent signals. 

However shortly afterward, Shekel pointed ont the ambignity problem in Ville’s work 
dehning the instantaneons parameters of a real signal [7, 8, 9]. As an example, consider 


x{t) = 3? 



t*;o(T)dT+(/)o] 


(7) 


There are an inhnite set of pairs ao(t) and a;o(^) for which x{t) may be eqnivalently de¬ 
scribed and hence an inhnite set of lA/IF parameterizations. Shekel pointed this ambignity 
ont “with the hope of banishing it [IF] forever from the dictionary of the commnnication 
engineer.” Others snch as Hnpert, snggested that despite the ambignity, the concept of in¬ 
stantaneons parametrization of real signals was still nsefnl and conld be applicable [10, 11]. 

In order to constrain the ambignity problem so that a nniqne complex extension can 
be jnstihed for a real signal, Vakman proposed three conditions tied to physical reality 
[12, 13, 14]: 1) amplitnde continnity, 2) phase independence of scaling and homogeneity, 
and 3) HC. With these conditions, Vakman showed that the nniqne complex extension is 
given by the rnle y{t) = 

More recently, other anthors have proposed other conditions to constrain the ambignity, 
snch as bonnded amplitnde and bonnded IF variation, leading to alternate lA/IF solntions 
[15, 16]. Finally, we point ont that a nnmber of other methods have also been introdnced, 
withont explicitly stating the conditions. Vakman has shown that all these methods violate 
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one or more of the conditions that he proposed and almost all retain HC [17, 18]. Like many 
authors, we believe Vakman’s conditions justifying Gabor’s use of the HT are reasonable 
and sensible for several special cases of z{t). On the contrary, we believe that Vakman’s 
conditions, in particular HC, can be too restrictive for more general z{t) and in the worst 
case, may lead to incorrect results and interpretations. 

Part I is organized as follows. In Section 2, we review the HT and its motivations and the 
AS. In Section 3, we give a proof showing we can still maintain analyticity without HC and 
also give a proof that no universal rule for the quadrature exists. In Section 4, we provide 
an example of a LSA problem and several solutions. Finally, in Section 5 we summarize 
Part I. 

2. The Hilbert Transform and Analytic Signal 

Although there exist several methods for estimating the instantaneous parameters, use 
of the HT dominates science and engineering. The HT of x{t) is given by 

CXO 

^ --f ^dT (8) 

ttJ x-t 

— OO 

where indicates the Cauchy principle value integral [19, 20]. The three main motivations 
for use of the HT are: 1) Vakman’s physical conditions, 2) analyticity of the resulting 
complex extension, and 3) ease of computation via Gabor’s QM. In this section, we review 
the motivations. 

2.1. Vakman’s Physical Conditions 

Vakman proposed conditions in order to constrain the ambiguity in choosing the complex 
extension [12, 13, 14, 21, 18]. 

Condition 1 Amplitude Continuity: Simply stated, amplitude continuity requires that 
the lA, pit) is a continuous function. This implies that the rule, y{t) = C{x{t)} must be 
continuous, i.e. 

C{x{t) + ew{t)} —)■ C{x{t)} for ||er(;(f)|| —0. (9) 

Condition 2 Phase Independence of Scaling and Homogeneity: Let x{t) have a complex 
extension, z{t) = p(f) exp[j0(f)]. Then for a real constant c > 0, cx{t) has associated 
complex extension Zi{t) = [cp(f)] exp[j0(f)], i.e. only the lA of the complex representation 
is affected and 0(f) and Q{t) remain unchanged. This implies that the rule for performing 
the complex extension is scalable 


C{cx{t)} = cC{x{t)}. (10) 

Conditions 1 and 2 force the operator to be linear [18]. 

Condition 3 Harmonic Correspondence: Let x{t) = Oq cos{uot + (fo), then HC forces the 
complex extension, 

z{t) = (11) 

4 



where we note the lA and IF must be constant. This implies that 

£{ao cos(a;ot + (/>o)} = Oo sin(a;ot + 0o) (12) 

and z{t) is a Simple Harmonic Component (SHC) with positive IF. 

As Vakman showed [12], the HT is the only operator that satishes the above conditions 
and as a result, the HT (or Gabor’s practical QM implementation) is viewed as the correct 
way to complex extend a signal. With his work, Vakman was able to refute most objections 
as to whether use of the HT is physically justihed. 

Condition 4 Phase Continuity: In addition to real signals having an ambiguity in instan¬ 
taneous parameterization, complex signals also have an ambiguity: although we construct 
z{t) using a rule for y{t), we want the lA/IF pair p(t), Q{t) for signal analysis, and there 
can be ambiguity in this coordinate transformation [22] . There exist at least two choices for 
resolving this ambiguity when it arises: 

• Condition 4a Positive lA: p{t) > 0 V h 

• Condition 4b Phase continuity: the phase 0(t) is a continuous function. 

Although condition 4a is the traditional choice, we advocate condition 4b to ensure the IF 
is well dehned. This is apparent in (3) where the lA may be negative. 

2.2. Cabor’s Quadrature Method 

The HT is an ideal operator that in practice is not physically realizable. Gabor’s QM 
provides a frequency domain approach that under certain conditions is equivalent to the 
HT. The steps in Gabor’s QM are: 

1. decompose x{t) into SHGs, i.e. compute the Fourier spectrum 

2. double the magnitude of the non-negative frequency components 

3. negate the negative frequency components. 

Gabor’s QM results in a complex signal formulated in terms of non-negative spectral fre¬ 
quencies, 

K-l 

m = E afcexpjj [wfcf+ 0fc]} (13) 

k=0 

i.e. we have decomposed the signal into SHGs each having constant lA, and (non-negative) 
constant IF, Uk- By comparing the expressions for z{t) in (13) and (1), we can effectively 
collapse K SHGs into a single AM-FM component and then obtain an lA/IF pair. Im¬ 
plementing Gabor’s QM using a FT is very convenient and is one reason for the method’s 
popularity. 
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2.3. Analyticity of the Analytic Signal 

Unfortunately, the word “analytic” has two distinct meanings when used in signal pro¬ 
cessing and in addition, another meaning in mathematics. A signal is said to be analytic if 
it consists only of non-negative frequency components [1, 23]. The analytic signal refers to 
the complex extension of a real signal using the HT, i.e. Gabor’s method [24, 25, 1, 26, 23]. 
If z(t), where t = ^ + jA, is an analytic function then the real and imaginary parts of the 
complex function, 

Z{i) =u{t,T)+jv{t,T) (14) 

satisfy the Cauchy-Riemann (CR) conditions 

d d d d 

—M(t,r) = —n(t,r) and —m (t, r) =-—n (t, r). (15) 

For the AS z{t) = x{t) + i'H{x{t)}, such that t ■(— f t -|- jr, the complex function 
i(t) = uit,T) + jv{t,T) is an analytic function [27]. Hence the reason for calling a HT- 
extended signal the AS [24, 1]. 

3. Relaxing the Condition of Harmonic Correspondence 
3.1. On Harmonic Correspondence 

We can understand Vakman’s motivation for tying HC to physical reality by considering 
the differential equation 

^z{t) + ujlz{t) = D (16) 

which describes many ideal systems, e.g. a LC circuit or mass/spring model. The solution 
to (16) is the SHC in (11). 

Any deviation from this ideal case, requires modihcation of the differential equation. 
For example, when the differential equation describes a circuit and resistance is included or 
describes motion and damping is included, the equation becomes 

d 2 / \ 

= 0 (17) 

where c is a constant. In this case, the solution includes an AM term and has the form 

z{t) = (18) 

which is not a SHC [3]. As another example, when the differential equation is further mod- 
ihed to include time-varying coefficients, non-linearities, or partial derivatives with respect 
to r or spacial variables, the resulting solution may include both AM and FM terms, which 
is also not a SHC [28, 18, 29]. 

While most authors believe HC is a reasonable condition to describe physical phenomena, 
we advocate that this condition is overly constraining. For many analysis problems, this can 
lead to incorrect interpretations because real physical systems always deviate from the ideal 
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case. By not assuming HC, we gain a degree of freedom in our analysis that allows us to 
construct other complex extensions that may be better suited to describing the underlying 
physical phenomena associated with the signal. We believe that any attempt to find a unique 
rule to infer z{t) of the form in (1) from x{t) in the form of (2) is fundamentally flawed and 
that no such universal rule can exist. 

3.2. On Analyticity of the Complex Extension 

We believe that the use of the term “analytic” has in most cases become too restrictive in 
signal processing. To wit, one may falsely believe that the HT is the only way to complex- 
extend a real signal in order to result in an analytic function, when time is considered 
complex. Other complex extensions of real signals can be constructed that result in analytic 
functions. Choosing y(t) = ^{{xit)} to obtain the AS z(t) ensures z(X) is an analytic 
function. However, there are other choices for y{t) such that z{i) is an analytic 

function. 

Theorem If we do not assume HC, there exists at least one choice for the quadrature, 
y(t) 7 ^ 'H{x(t)} that results in z(t) being an analytic function. 

Proof. Let x{t) = aocos{uot) and choose the quadrature such that 

z{t) = aocos{uot) + jaaosm^fuot) (19) 

with real a and f. The complex function is given by 

^(t) = aocos(a;o[t + jr]) + jaaosm{/3uo[t + jr]) (20) 

where 

u{t, t) = oo cos(a;ot) cosh(a;oT) — aa^ cos{(3uot) sinh(/9a;o'7‘) (21) 

and 

v{t, t) = aao sm(fujot) cosh(/3a;or) — oq sin(a;ot) sinh(a;oT). (22) 

It can easily be shown that this choice leads to z(t) satisfying the CR conditions and hence 
is an analytic function. Any choice oi a f 1 oi f 1 does not imply HC. □ 

Although Gabor’s QM provides a simple rule to obtain an lA/IF pair that parameterizes 
x{f), it may not address the problem of obtaining the latent signal from the observation. 
The heart of the problem is that no single rule can determine the latent signal from the 
observation for every possible latent signal because more than one complex signal map to 
the same real signal under the real operator. Although a unique rule can be constructed to 
work for a particular z{f), the same rule cannot generally be used. 

Corollary No unique rule for the quadrature, y{f) = £{x(t)} exists to obtain the latent 
signal, z{t) from the observation, x{f) = 5R{2:(t)} for all z{t). 

Proof. Consider the latent signal, z{t) of the form in (19). Assume a unique rule, y{t) = 
C{x{t)} = oo^osin(/3oa;o^) exists to obtain z{t). This implies a unique z{t) = aocos{uot) + 
jaoOo sin(/3oa;o^)- li ao a or fo (3 then z{t) z{f) and /!{■} is not unique. □ 
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3.3. On Harmonic Conjugate Functions 

If z(t) = u(t,T) + jv(t,T) is an analytic function, then u(t,T) and v(t,T) are unique and 
are harmonic conjugates [27]. Even though x(t) = u(t,0) and y(t) = v(t,0) in (1) this does 
not imply that u(t,T) = a;(t) and v(t,T) = y(X), but rather u(t,T) and v(t,T) each contain 
terms from both a;(t) and |/(t): 

^(t) = x{i)+jy{i) 

= [xR{t, r) + jxi{t, r)] + j [ynit, r) + 3 yi{t, r)] 

= [xR{t, r) - yi{t, r)] + j [xi{t, r) + ynit, r)] 

= u{t,T) + jv{t,T) (23) 

where R, I denote the real and imaginary parts, respectively. This is easily seen in (21), 
where a and (3 are present in u{t, r) despite a and [3 appearing only in y{t) in (19). Thus the 
problem of Ending y{t) cannot be solved by finding a harmonic conjugate of u{t, r) because 
y{t) must be known to obtain u{t, r). 

3 . 4 . On AM-FM Demodulation 

The HT is widely used as a demodulation algorithm for AM-FM signals. This is typically 
justified as a valid approach because of Vakman and also Bedrosian’s theorem.^ The HT 
can be used to determine the lA/IF for a small subset of latent signals, i.e. those with 
HC. The HT can also be used to closely approximate the lA/IF for latent signals whenever 
R{x{t)} ^ y(t) [8]. However, the HT cannot be used in general to obtain the lA/IF of all 
latent signals. 


4. Example 


Recall the LSA problem: given the observation x{t) = 3fJ{z(f)}, find z{t) or more strictly 
the instantaneous parameters, p{t) and fl{t). In this section, we will give three solutions to 
a LSA problem where two of these solutions have y{t) ^ 'R{a;(t)}, as illustrated in Fig. 1(b). 

Suppose we have three ideal systems: a frequency modulator, an amplitude modula¬ 
tor, and a Linear Time-Invariant (LTI) system in the steady state. We observe a triangle 
waveform x{t) that could have come from any of the three ideal systems. What is the 
corresponding latent signal z{t) or more strictly the instantaneous parameters, p{t) and 
D{t)l 

Let the observation be the periodic, even-symmetric triangle waveform where one period 
is given by 


f —2AijjQt/'K + A, 0<f<T/2 
I 2AuQt/TT + A., —T/2<f<0 


(24) 


with amplitude A, period T, and fundamental frequency wq, as illustrated in Fig. 2. 


^If the product consists of a low-frequency factor, l{t) and a high-frequency factor, h{t) of non¬ 
overlapping spectra, then 'H{l{t)h{t)} = [30, 31, 18]. 
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Figure 2: One period of the triangle waveform, x(t) in (24) with amplitude A and ujq = 1. 


4-J- Solution Assuming Harmonic Correspondence 

If we assume HC, i.e. the ideal system is LTI in the steady state then 

Zi{t) = x{t) + j'H{x{t)} (25a) 

= ^2(2fc + l)2 ^ Y. ^2(2fc + i)2 [(2^ + 1)^0^] (25b) 

where the lA is 

8/1 

pit) = —doit) (26) 

TT^ 

and the IF is 

Clit) =020 + —Moit) (27) 

where ao(t) and Moit) are related through 

OO ^ 

4-2. Solution Assuming Constant IF 

If we assume constant IF flit) = uo, i.e. the ideal system is an amplitude modulator then 

Z2it) = pit) cosiuot) + jpit) sin(i:not) (29) 

with the lA given by 

pit) = xit)/cosiuiot). (30) 

In this solution, y(t) = pit) sin(a;ot) ^ 'H{/3(f) cos(a;ot)}, because Bedrosian’s theorem cannot 

be applied. 
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4-3. Solution Assuming Constant I A 

If we assume constant lA p{t) = A, i.e. the ideal system is a frequency modulator then 
Z 3 {t) = A cos [uot + Mo{t)] + jAsin [uot + Mo{t)] (31) 

with IF given by 

Cl{t) = cuo + —Mo(t) (32) 

where 

Mo(t) = arccos [a;(t)/A] — Wot. (33) 

As in the solution for constant IF, y{t) = A sin [u^t + Mo(t)] ^ "HjAcos [u^t + Mo(t)]}. 

5. Summary 

We have presented the LSA problem and reviewed the motivation for the use of the HT 
and AS in signal analysis. We proved the HC condition is not necessary for analyticity and 
that no unique rule for the complex extension exists to obtain the latent signal. It was 
argued that by relaxing the HC condition, which forces the use of the HT, we gain alternate 
choices for the latent signal. Although the HT is widely used for demodulation of AM-FM 
signals, it cannot be used to obtain the lA/IF of all latent signals. In a strict sense, Gabor’s 
QM can only be used when the latent signal is a superposition of SHCs with non-negative IF 
and can approximate the latent signal whenever ^/{^(t)} ~ y{t), e.g. communication signals 
where Bedrosian’s theorem is approximately satished. 
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Part II. Hilbert Spectral Analysis Using Latent AM—FM Components 

[W]ell-known methods ... can no longer be applied when the freqnency itself is 
made a fnnction of the time. Even the concept of “instantaneons-freqnency”... 
is of a somewhat arbitrary bnt nevertheless highly usefnl natnre, is often mis- 
nnderstood and even misinterpreted .... The only way at present available to 
solve these and similar problems is to go back to the very first and fnndamental 
principles. This implies a theoretical treatment beginning with the differential 
eqnations of the problem concerned. Unfortnnately these eqnations can seldom 
be solved in terms of well-known fnnctions, snch as real or complex exponentials 
[simple harmonic components]. I think it is precisely to this fact that the main 
difficnlties can be traced.—Van der Pol [5] 

In the second part of the paper, the central theme is Hilbert Spectral Analysis (HSA) 
using a superposition of latent AM-FM components. We present HSA as a generalized 
LSA problem. In the general problem, we seek a representation of z{t) consisting of a 
superposition of latent components, i.e. a multicomponent model. Furthermore, rather than 
seeking a single lA/IF pair for z{t), we seek a set of lA/IF pairs each associated with the 
components. Although time-frequency analysis has been extensively studied the use of a 
generalized AM-FM model for this analysis, without the HC condition, has never been 
proposed. Using this model leads to non-unique signal decompositions due to the theorem 
proved in Part I. However, by imposing assumptions on the form of the AM-FM component, 
a unique parameterization in terms of lA and IF can be obtained. This analysis requires 
abandoning Gabor’s complex extension as advocated in Part I and instead allowing the 
assumptions to imply the complex extension. This model enables us to analyze signals 
with very few restrictions resulting in alternate and possibly more useful decompositions, 
especially for nonstationary signals. In this part, AM-FM modeling and HSA theory is 
presented. Examples using HSA are given and a visualization of the Hilbert spectrum is 
proposed. 

1. Introduction 

It has long been known that simple harmonic analysis of nonstationary signals may lead 
to incorrect interpretations of an underlying signal model [9]. In Gabor’s seminal paper, 
“Theory of Gommunications” he writes: 

The greatest part of the theory of communication has been built up on the 
basis of Fourier’s reciprocal integral relations .... Though mathematically this 
theorem is beyond reproach ... even experts could not at times conceal an uneasy 
feeling when it came to the physical interpretation of results obtained by the 
Fourier method. After having for the hrst time obtained the spectrum of a 
frequency-modulated sine wave, Garson wrote: ‘The foregoing solutions, though 
unquestionably mathematically correct, are somewhat difficult to reconcile with 
our physical intuitions ....’ [6] 
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As Gabor noted, Carson in 1922 was the first to rigorously study a FM signal and 
realize that the frequency components obtained from such a nonstationary signal, do not 
describe the physical system in a meaningful way [32], A similar argument regarding an 
AM signal was made by Priestly, who wrote, “...a nonstationary process in general cannot 
be represented in a meaningful way by the simple Fourier expansion” [33]. His example, 

x{t) = cos{uot + 0o) (34) 

and its Fourier Transform (FT) consists of two Gaussian functions centered at frequencies 
icuo- Thus, the FT contains an infinite number of SHCs. This interpretation of the underly¬ 
ing signal model may be incorrect. For example, an alternative signal model consists of just 
two components at constant frequencies ic^o, with each component having a time-varying 
amplitude, (A/2)e“* . Mathematically, these two representations are equally valid and 

correspond to different families of basic functions used for representation [9]. 

For nonstationary signals, SHCs lose effectiveness and thus the idea of time-varying 
components^, i.e. components with time-varying lA and IF has arisen in order to account 
for the non-stationarity [33, 9]. However, the concept of IF is not without its own controversy 
primarily due to four reasons: 

1. There is an apparent paradox in associating the words “instantaneous” and “fre¬ 
quency” because frequency usually defines the number of cycles undergone during 
one unit of time [9]. 

2. Without assumptions, instantaneous parameterizations of a signal are not unique [7, 
9, 1, 16]. 

3. The commonly accepted definition of IF, i.e. phase derivative of Gabor’s AS is correct 
for only a limited class of signals [1] . 

4. Although different quantities, harmonic frequency and IF are often confused likely due 
to the term “frequency” attached to both [34, 35]. Harmonic frequency is a special 
case of IF and is only equivalent under the assumption of SHGs. 

Several authors over the previous decades have shown that problems and paradoxes exist that 
are related to the definition of IF [7, 34, 16, 1, 35, 36, 21, 33, 37, 38, 39, 40, 41]. However on 
the whole, these problems seem to have been forgotten or ignored. In 1937, Garson and Fry 
formally defined the IF based on the phase derivative of a complex FM signal assuming the 
FM message is slowly-varying [4] . This assumption, while perfectly valid in communication 
theory, implies a narrowband component when used in signal analysis. 

In order to exploit the mathematical convenience of using complex exponentials in signal 
analysis, Gabor introduced a QM for obtaining a complex extension of a real signal [6]. Ga¬ 
bor’s QM assumed positive IF SHGs and was shown to be equivalent to the HT [6]. Although 


^Some authors, such as Ville [24], refer to time-varying components as “instantaneous spectra” or more 
generally as functions of time that give the structure of a signal at a given instant. 
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Gabor’s QM provides a unique method for complex extension, as Cohen pointed out without 
strict adherence to the assumption, this results in many counter-intuitive consequences [1]. 
Ville defined the IF of a real signal by using Gabor’s complex extension and then Carson’s 
definition of IF [9]. By defining the IF as the derivative of the phase of Gabor’s AS, Ville 
was able to show the average harmonic frequency is equal to the time average of the IF 
[40]. He then formulated the Wigner-Ville Distribution (WVD) and showed that the first 
moment of the WVD with respect to frequency yields the IF. Using Gabor’s AS to extend 
a real signal results in a number of very convenient relations. As a result, Gabor’s QM is 
almost universally viewed as the correct way to dehne the complex signal and subsequently 
the correct way to dehne lA, IF, and phase for real signals despite the inherent assumption 
of HC [1]. 

In our research, we have found that the dogmatic use of Gabor’s complex extension does 
not provide the necessary hexibility for modeling nonstationary signals. For many problems, 
as we will demonstrate, it can be advantageous to allow the assumptions of an underlying 
signal model to imply the complex extension. As a result, assumptions must be made on a 
per problem basis in order to force a unique decomposition and hence, in order to properly 
estimate the latent signal and its components. This leads to a decomposition into complex 
AM-FM components rather than SHCs. As a result of using AM-FM components, we 
revert back to Carson’s original dehnition of IF and by further relaxing HC and embracing 
the resulting non-uniqueness, the associated problems and paradoxes related to IF can all be 
resolved. 

Part H of this paper is organized as follows. In Section 2, we propose the AM-FM 
model. In Section 3, we propose a general Hilbert spectrum based on the AM-FM model 
and examine several familiar specializations. In Section 4, a frequency domain view of LSA 
is presented. In Section 5, we discuss subtleties of complex extension of real signals assuming 
AM-FM components with relaxed HC. In Section 6, we provide examples of HSA and in 
Section 7 propose a visualization of the Hilbert spectrum to aid in the interpretation. In 
Section 8, we summarize Part H. 

2. The AM—FM Model 

2.1. Definitions and Assumptions 

Our goal in HSA is to decompose a signal into complex AM-FM components, each of 
which yield lA and IF estimates that are matched to some criteria. That criteria could be 
some physical or perceptual observation or simply some intuitive notion of an underlying 
signal model. The decomposition problem is illustrated in Fig. 3 where we have extended 
the LSA problem shown in Fig. 1 to show models of the latent signal that in general, are 
composed of a set of latent AM-FM components. 

We propose to define the AM-FM model for z(t) as a superposition of K (possibly 
inhnite) AM-FM components 

K-l 

m - E ^|Jk{t;ak{t),Uk{t),(j)k) (35) 

k=0 
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Harmonic 

Correspondence 


Figure 3: A set diagram for the HSA problem. In the LSA problem, many latent signals z{t) map to the 
same observed signal x{t) under the real operator. In the HSA problem, many component sets {tjjk{t)}, k = 
0,1,... ,K map to the same latent signal z{t) under superposition. The goal in HSA is to properly choose 
{i’kit)} given x{t). 


where the AM-FM component is dehned as 




akit)exp < j 


Uk{x)dx + (j)k 


Sk{t) +jcrk{t) 


(36a) 

(36b) 

(36c) 


parameterized by the lA afc(t), IF 0 Jk(t), and phase reference (pk- We assume the observed 
real signal x{t) is related to the latent signal z{t) by (2). 

When convenient, we will drop the k denoting the parameters of the kth component 
for notational simplihcation and unspecihed references to lA, IF, and phase refer to ak{t), 
oJkit), and 6 kit) rather than pit), fl(t), and 0(t), respectively. We note that most of the 
paradoxes related to IF are a result of not correctly interpreting between these variables 
in the multicomponent case, because in the monocomponent case they are equivalent. The 
geometric interpretation of the AM-FM component in (36) is illustrated with the Argand 
diagram in Fig. 4. The AM-FM component can be visually interpreted as a single rotating 
vector in the complex plane with time-varying length and time-varying angular velocity. 

Carson expressed the phase 6*(t) in terms of a carrier frequency ujc and FM message m(t) 
as [4] 

t 

Oit) = Uct + y mix)dx (37) 

0 

assuming the modulation index, X < Uc and |m(f)| < 1. We choose to parameterize the 
phase as 


9 it) = uJct + 


mir)dx + (p 


— OO 


(38) 
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Figure 4: (a) Argand diagram of an AM-FM component in (36) at some time instant. Each component, 
( — ) is interpreted as a vector: the lA a{t) is interpreted as the component’s length, the phase 0{t) 
is interpreted as a component’s angular position. Although not shown, the IF uj{t) is interpreted as a 
component’s angular velocity and phase reference (j) is interpreted as an initial condition. The real part of 
the component s{t) ( — ) and the component’s quadrature a{t) ( — ) are interpreted as orthogonal projections 
of We have included an example path ( — ) taken by tpit). (b) Argand diagram of the latent signal z{t) 
in (1) at some time instant, composed of a superposition of components shown in red. The latent signal 
( — ) is interpreted as a vector: the lA p{t) is interpreted as the latent signal’s length, the phase 0(f) is 
interpreted as a latent signal’s angular position, and although not shown the IF 17 (f) is interpreted as a 
latent signal’s angular velocity. The real part of the signal x{t) ( — ) and the signal’s quadrature y{t) ( ) 

are interpreted as orthogonal projections of z{t). We have included an example path ( — ) taken by z(t). 

which is more general than (37), avoids an arbitrary normalization on m(t), and does not 
impose an upper bound on the IF. The phase can also be expressed in terms of Uc and Phase 
Modulation (PM) message M{t) as 

9(t) = ojct + M{t) + 0 (39) 

where the relationship between the messages is given by 

t 

M{t) = J m{x)dT. (40) 

— OO 

In the AM-FM model, we assume non-negative IF 

uj(t) = ^Oit) = ujc + mit) > 0 V t. (41) 

at 

In Section 5, we discuss situations in which this assumption may be relaxed. Finally, we 
assume without loss of generality, that the phase and carrier references are taken at t = 0, 
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i.e. 

0 

J m{t)dt = M{0) = 0 (j) = 9(0) and Wc = ti;(0) — m(0). (42) 

— OO 

2.2. Monocomponents and Narrowband Components 

The concept of a signal being composed of one (mono-) or more (multi-) components 
and whether those components are narrowband or wideband has caused much confusion 
in existing literature. There is no clear agreed upon dehnition for a monocomponent [42], 
however, a few dehnitions have been proposed [43, 1, 23]. Some authors describe a mono¬ 
component as dehned by a single ‘ridge’ in time and harmonic frequency, corresponding to 
an elongated region of energy concentration [43, 44, 1]. Cohen agrees that if a component is 
well localized, the crest of the ridge corresponds to the IF and further states that the width 
of the ridge depends on the energy spread, or instantaneous bandwidth of the component 
[44, 1]. A multicomponent signal may then be dehned as any signal which is the sum of two 
or more monocomponents which can only occur if the separation between the ridges is large 
in comparison to bandwidth of the components [43, 44, 1, 45]. It has also been noted that 
a signal may be monocomponent at some time instances and multicomponent at other time 
instances [44]. 

With these descriptions of a monocomponent, we point out that the dehnitions and 
concepts of lA and IF of a complex AM-FM component as dehned by (36) are by no 
means justihed only for narrowband components. Also, as pointed out by Cohen [1], a 
multicomponent signal is not dehned by the ability to express a signal as a sum of parts, 
because there are an inhnite number of ways to write a signal as a sum of parts. Rather, it 
is the nature of the parts in relation to themselves and the signal which determines whether 
the decomposition is of interest [46]. We strongly believe that the previous description of a 
monocomponent is unnecessarily restrictive, limiting the monocomponent to be narrowband. 
Clearly, the lA and IF in (36) are well-dehned for wideband components. 

Simply put, a perfectly-valid component can be dehned with (36) while foregoing the 
restrictive narrowband dehnitions of Carson [4], Ville [24], Boashash [9, 43], and Cohen [1]. 
Such wideband components may result in multiple ridges which may be broken down into 
narrowband components. However, for many problems a wideband component, such as the 
AM-FM component in (36), can be much more appropriate than multiple, narrowband com¬ 
ponents. Further, the appearance of structure in the Fourier spectrum can be viewed as an 
indication that a wideband component is present in the signal. For narrowband components 
the Fourier and Hilbert spectra can be quite similar, but for wideband components these 
can be quite diherent. 

3. Hilbert Spectral Analysis 

Huang’s original dehnition of the Hilbert spectrum uses Empirical Mode Decomposition 
(EMD) to determine a set of Intrinsic Mode Functions (IMFs) which are individually de¬ 
modulated with the HT to obtain {ak(t)} and {uk(t)} [42]. This analysis is also known as 
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the Hilbert-Huang Transform (HHT) [42], This dehnition can be generalized, by recognizing 
that IMFs are a class of AM-FM components and that other decomposition and demodn- 
lation methods exist for obtaining the instantaneons parameters, {ak{t)} and This 

generalization can lead to a more powerfnl and usefnl signal analysis techniqne as we will 
demonstrate. 

Towards a more general definition, we define the Hilbert spectrnm as a representation or 
characterization of a signal by an instantaneous spectrum which is parameterized by AM-FM 
components as in (36). With this definition the Hilbert spectrnm is not nniqne, therefore 
two problems mnst be addressed in order to obtain a nniqne solntion: 

PI: Determine the appropriate complex extension to obtain an AM-FM decomposition 
with meaningfnl interpretation. 

P2: Estimate {ak{t)} and {uk{t)}, 0 < k < K — 1. 

The instantaneons parameters are defined on a complex signal, therefore to perform 
analysis on a real signal it mnst be extended to the complex space. There are cases where 
the extension is straightforward, i.e. simple harmonic motion, or in commnnications where 
assnmptions can be matched at the transmitter and receiver, or other cases where the 
assnmptions are implicit, i.e. Fonrier analysis. However, in the general case, signal decom¬ 
position and component demodnlation are ambignons and dne to this, assnmptions mnst 
always be made in order to arrive at a nniqne solntion. We will address snbtleties of complex 
extension of real signals assnming AM-FM components with relaxed HC in Section 5. 

The class of the solntion obtained is dependent on the natnre of the assnmptions. Fig. 5 
illustrates the various classes of components under particular assumptions including the 
SHC, AM, FM, and AM-FM components as well as the IMF (discussed further in Part 
HI). The AM-FM component in (36) may be considered a generalization of all of these 
components. In the next subsections, we illustrate several forms of the AM-FM component 
and AM-FM model under various assumptions. 


AM-FM Component 



Intrinsic Mode Function 


AM Component 




Harmonic Component 

FM Component 







Figure 5: The AM-FM component in (36) may be considered a generalization of other well-known compo¬ 
nents. The familiar harmonic component may be considered special cases of the AM component and FM 
component. Huang’s IMF is a special case of the AM-FM component. 
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3.1. Two Conventions for Relating the Real Observation to the Latent Signal 

As noted by Gabor in [6], there are two conventional ways to relate a real signal x{t) to 
a complex signal z{t). The hrst convention is 

x{t) = 

= ^{x{t) + jy{t)} 

which can be thought of as a single vector as in Fig. 4 and the second convention is 

x(t) = [z(t) + z*(t)] /2 (44a) 

= [x{t) + jy{t) + x(t) - jy(t)]/2 (44b) 

which can be thought of as two vectors. If z(t) is a superposition of AM-FM components, 
the hrst convention can be viewed as a single rotating vector per component, while the 
second convention can be viewed as a pair of vectors per component rotating in opposite 
directions. 

Equation (44) is implicit in standard Fourier analysis, with the interpretation that the 
latent signal is real-valued and its spectrum^ consists of Hermitian symmetric, positive and 
negative frequency components. The convention in (43) is adopted for the AM-FM model in 
(35), with the interpretation that only the real part of the signal, x{t) is observed. Inherent 
in both these conventions is the ambiguity associated with the quadrature y{t), i.e. an inhnite 
number of choices exist for y{t) that lead to the same x{t). Thus, we view the quadrature 
y{t) as a free parameter that can be chosen to achieve alternate interpretations. 

We use the term quadrature in a general sense. For mathematical and physical reasons, 
the quadrature signal is most often chosen offset in phase by one-quarter cycle /2 radians) 
relative to the real signal. More specihcally, if we assume a single SHC, z{t) = fjoit] ag, ovg, 0o) 
then the AM-FM model is given by (11). Although choosing the quadrature signal given 
the real signal is trivial for the SHC, the same is not true for the component in (36) where 
the assumptions ao(t) = ag, and LVg(t) = ug are removed, i.e. the concept of a quadrature 
signal for the AM-FM component is ambiguous—as will be demonstrated in the following 
sections. 

3.3. Simple Harmonic Component 

The simplest form of the AM-FM model has K = 1, ag{t) = ag, and u:g{t) = ujg in which 
case the AM-FM component in (36) becomes 

V’o(i; ao,a;o, (fg) = (45) 

We recognize the signihcance of this form as representing simple harmonic motion as we 
discussed in Part I as Vakman’s Condition 3. 


^In this work, use of the unspecified word “spectrum” strictly means Fourier spectrum and not “Hilbert 
spectrum.” 


(43a) 

(43b) 


18 



3.3. Superposition of Simple Harmonic Components 

Perhaps the most familiar form of the AM-FM model has K = oo, ak{t) = Ok, and 
LJkit) = kujQ 

OO 

z{t) = (46) 

k=0 

and has real observation given by (2) repeated here as 

x{t) = 5R 

OO 

= ^ afc cos{kuot + (fk) 

k=0 

oo 

= Aq + Ei^ k cos{kuot) + Bk sin(/ca;ot)] 

k=l 

Equation (47) is recognized as a Fourier Series (FS) with a particular convention used for 
the complex extension [47, 48]. Note that (fk is required for proper synthesis since without 
it, Ak = Ok and Bk = 0 and the resulting x{t) is always even, which may not be true. The 
FS can be considered the simplest of AM-FM superposition models corresponding to the 
assumption of SHCs. 

The FT is the limiting form of the FS as the fundamental, Wq —t 0 [49]. When viewing 
the FT as a special case of the AM-FM model, subtle and important observations can be 
made: 

• In the decomposition problem PI, the AM-FM model corresponding to a separable 
solution with X{uj) (not time-varying) and (constant frequency), is the FT, i.e. the 
inverse FT is a superposition (inner product) of X{u) and 

• The FT yields a solution that can be described as a superposition of non-time-varying 
components described by constant state parameters, ak(t) = Ok, oJkit) = kuo, and (fk 
which corresponds to a constant model state. However, this does not in any practical 
way constrain the real superposition x{t). Although analysis of a time-varying system 
can be accomplished using a constant model, like the FT, in most cases this forces 
K = oo which may not accurately describe a physical system. 

• HSA can be considered a generalization of Fourier analysis [42]. However, rather than 
a generalization of Fourier analysis obtained by generalizing a kernel function (Gabor 
transform, Wigner-Ville distribution, wavelet transform, etc.) [50, 51, 52, 53, 54, 1, 55], 
an AM-FM component can be viewed as allowing for additional degrees of freedom 
in the bases of Fourier analysis. Although it may be convenient to think of (36) as a 
basis for HSA because they span the entire Hilbert space, they do not meet the linear 
independence property of a formal basis [56]. 


(47a) 

(47b) 

(47c) 
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3.4- The AM Component 

One way to generalize the SHC in (45) is to relax the constraint of constant amplitude 
which leads to an AM component 

(48) 

The simple AM component was originally developed in the context of communication theory 
[57, 58, 59], 

3.5. Superposition of AM Components 

Another form of the AM-FM model is a superposition of AM components that have 
frequencies as integer multiples of a fundamental 

OO 

(49) 

k=0 

with real observation given by (2). 

The AM superposition was first conceived by Gabor and lead to the development of the 
Short-Time Fourier Transform (STFT) [6] 


Xt{uj) = f x{x)h{T: — t)e 

(50a) 

J 

^ 1 f 


Xi^{t) = — / X{cy)H{u — (Cje^^^dci. 

2vr J 

(50b) 


where h{t) ^ H{u) is a window function. There are two possible interpretations of the 
STFT in context of the AM-FM model; 

1. The window function h{t) is applied to the signal x{t) in order to use the FT, i.e. simple 
harmonic analysis at each t. This interpretation couples H{u) to X{u). 

2. The window function h{t) is applied to the and is effectively an AM superposition 
analysis where the window is an assumption on ak{t) in the AM-FM model. This 
interpretation couples H{u) to 

The second interpretation provides an important and alternate view of the STFT. What 
is typically viewed as imprecision in the ability to compute time-/reg'uenc^ parameters in 
the first interpretation can alternatively be viewed in the second interpretation as a precise 
ability to compute time-component parameters using the modified basis, h{t — T)e“-^‘^*. 

3.6. FM Component 

Another way to generalize the SHC in (45) is to relax the constraint of constant frequency 
which leads to a simple FM component 

t 

j[ f uJoOdr+rpo] 

'ilJo{t;ao,Uo{t),(f)o) = aoc 
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(51) 


The FM component was also developed in the context of communication theory [32, 4], The 
FM component has been used in signal synthesis where Chowning observed that a single, 
wideband FM component is perceived by the ear as spectrally rich, i.e. multiple SHCs [60]. 
This FM-synthesis method has been used in commercial audio synthesizers [61, 62], 

3.1. Superposition of FM Components 

Signal analysis using a superposition of FM components is not usually considered due to 
the loss of linear independence of the components in the model and the restriction of constant 
amplitude. For this reason, signal analysis using a superposition of AM-FM components is 
more useful than a superposition of FM components. 

3.8. Other AM-FM Models 

Alternatives to the proposed AM-FM model in Section 2 have been considered. However, 
common to these are restrictive assumptions which limit the utility as we highlight below. 
Previous AM-FM models for signal analysis/synthesis usually fall into one of three main 
groups: 1) HT [63, 64, 65, 66], 2) peak tracking/sinusoidal modeling [67, 68, 69, 70], and 3) 
Teager energy operator [71, 72, 73, 74, 75, 76]. However, some models exist that do not fall 
into any of these groups [77, 73]. A historical summary of AM-FM modeling is presented by 
Gianfelici [78]. A review of algorithms for estimating IF is presented by Boashash [79, 80]. 

The HT permits the unambiguous dehnition of lA, IF, and phase of any real signal 
(random or deterministic) [21]. Thus, it provides a direct means of performing AM-FM 
analysis, however, it requires that the quadrature signal be dehned in a consistent manner 
in all cases. Implicit in the HT is HC, as a result the more likely this assumption is true, 
the better the Hilbert-transformed signal approximates the true quadrature signal and the 
more likely the HT-based solution provides an accurate model of the underlying physical 
synthesis model [43]. Direct application of the HT does not however address the signal 
decomposition problem, as a result, methods for using the HT for decomposition have been 
proposed [65, 81, 66]. 

In peak tracking, one accepts the narrowband dehnition of a component and as a result, 
each component appears in the time-frequency plane as a single ridge of energy concen¬ 
tration. Thus, a signal can be parameterized by tracking its ridges in location, intensity, 
and possibly bandwidth. The time-frequency distribution is usually derived from a STFT 
[67] but a generalized time-frequency distribution can also be used [68]. Regardless of the 
time-frequency distribution used, the narrowband assumption is inherent. Extensions of 
the sinusoidal model have been proposed, such as the harmonic plus noise, and adaptive 
quasi-harmonic model [69]. 

The Teager Energy Separation Algorithm (ESA) provides a method of estimating the 
lA and IF of an AM-FM component, assuming HC. The Teager Energy Operator (TEO) is 
dehned as [82] 

F {x{t)} = [x{t)]^ — x{t)x{t) (52) 
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where x{t) and x{t) are the hrst and second time derivatives of x{t), respectively. The TEO 
applied to a single narrowband AM-FM component in (36), resnlts in [82] 

~ [ao(t)a;o(t)]^ (53) 


Assuming the modulating signals ao{t) and mo(t) are bandlimited, it can be shown that 


and 


ao{t) ^ 





Uo{t) ^ 






(54) 


(55) 


thus providing a method to estimate narrowband monocomponent lA/IF parameters [82, 
83, 84, 85], 

The AM-FM model in these groups all rely on a rigid narrowband component. Surpris¬ 
ingly, the use of wideband components is a well-known means of synthesis [86, 62, 60, 87]. 
However, analysis using the wideband components in (36) in the general form, has not been 
considered. It is our belief that the true power of HSA can only be fully recognized with 
the use of wideband components without the inherent assumption of HC. 


4. Frequency Domain View of Latent Signal Analysis 

In Part I, we introduced the LSA problem in the time domain. The frequency domain 
view of the LSA problem is to estimate the latent spectrum Z{u) from X(a;). Because x{t) = 
3?{2;(f)}, this imposes structure on the spectrum of x{t), namely X{u) = [Z{u) + Z*{—u)]/2^ 
i.e. Hermitian symmetry. As discussed in Part I, the conventional way to estimate z{t) from 
x{t) is via the HT as in (6). Equivalently in the frequency domain, the conventional way to 
estimate Z{u:) from X{u) is with Gabor’s QM. Using Gabor’s QM effectively forces the FT 
of (5) Y{u;) = —jsgn{u)X{oj), and Z{u) is always spectrally one-sided [70]. In the frequency 
domain, by relaxing HG we gain the freedom to choose Y(u) appropriately. 

As an example of this frequency domain view, consider the observed spectrum 

X (ca) = ao7i[6{u + Wq) + (5(ta — Wq)] (56) 

where (5(-) is the Dirac delta function. If we assume HG, 

Y{u) = -jaoTr[-6{u + Uo) + 6{u - Wq)] (57) 

and the corresponding latent spectrum is given by 

Z{cj)=X{u)+jY{u) 

= 2ao7i6{u — coq). 
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(58) 







If we do not assume HC, there are many choices for Yioj) leading to many possible latent 
spectra including, 

Z{(jj) = 2aovr(5(a; + wq), (59) 

Z{(jj) = aovr[(l - a)6{u + Wq) + (1 + a)5{u - Wq)], (60) 

and 

Z{u;) = ao7r[(5(a; + wo) + ^(ti^ — cuo) + aS{oj — (3u;o) — C(6{u; + l3uo)] (61) 

where a,/9 G M and (61) corresponds to the FT of (19). Because of the structure imposed 
on the spectrum X{u) = [Z{u) + Z*(—a;)]/2, the latent spectra in (58)-(61) all yield the 
same X{u;). 

In Table 1 we summarize the spectral structure of a complex signal depending on whether 
or not HC is assumed. If z{t) takes on only real values, then the spectrum has Hermitian 
symmetry (first row). If z{t) takes on complex values and has HC, then the Fourier spectrum 
is single-sided (second row). If z{t) takes on complex values and does not have HC, the 
spectrum does not have Hermitian symmetry. 

Table 1: Structure of the Fourier spectrum under the assumption of a model composed of Simple Harmonic 
Components (SHCs). For a complex signal that takes on only real values, the Fourier spectrum is Hermitian 
symmetric. For a complex signal with HC, the Fourier spectrum is one-sided.® For a complex signal without 
HC, the Fourier spectrum is not Hermitian symmetric. 


Class of z{t) FT Structure 


z(t) G M 

z{t) G C w/ HC 

z{t) G C w/o HC 


Z{uj) = Z*{-u;) 
Z{uj) =0, o; < 0 
zluj) ^ Z*{-u:) 


In the frequency domain view of LSA, Gabor’s QM can be interpreted as simply a 
method to distinguish Z{u:) from Z*{—u:) which works only if the harmonic frequencies 
of Z{u:) are all non-negative. This is a convenient “trick” and a generalization of this 
concept for the Hilbert spectrum is possible. As an illustration, consider Fig. 6(a) where the 
latent spectrum Z{oj) = Z 2 {uj) -|- Zi(uj) and the real signal under analysis has a spectrum 
[ZJ'(—cn) -|- Z 2 (—u) -f Z 2 (u>) + 2'i(a;)]/2. Applying Gabor’s QM, the complex signal clearly 
has the spectrum Z(lj) = Z 2 {oo) + Zi(uj). Next consider Fig. 6(b) where the latent spectrum 
Z(lj) = Z 2 (lj) + Zi(uj) and the real signal under analysis has a spectrum [Zl{—u) + Z 2 (lj) + 
Z^i^—uj) + Zi(a;)]/2. Applying Gabor’s QM, the complex signal has the spectrum Z^^—u) -f 
Zi{u) which is the incorrect spectral grouping. Finally, consider Fig. 6(c) where the latent 
spectrum Z{u)) = Z^^u) + Z 2 (lj) + Zi{u) and the real signal under analysis has a spectrum 
[ZJ'(—cn) -f -|- Z^l:^ij)--^rX 2 {LJ) + Zi(c(;)]/2; cancellation is due to symmetry 

of the latent spectrum. Applying Gabor’s QM, the complex signal has the spectrum Zi(u) 


®For convenience, we omit the case that Z(uj) = 0, a; > 0. 
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and not only has Gabor’s QM failed to yield the latent spectrnm bnt terms are missing 
entirely. 


5. Subtleties of the Hilbert Spectrum 

In the previons section, we considered the LSA problem where we assnmed SHCs bnt 
relaxed the HC assnmption. By relaxing HC, we also relax the non-negative IF constraint 
associated with it. In this section, we discnss the implications of assnming non-negative IF 
and relaxing the assnmption of SHCs. 

First, consider making the assnmption of non-negative IF and a model composed of 
SHCs. Table 1, second row shows this implies HC on z{t). Next, consider non-negative IF 
and a model not composed of SHCs. Contrary to Table 1 rows two and three where SHCs 
are assnmed, there can exist models with non-negative IF bnt withont HC. We illnstrate 
the existence of snch models with the following example. Consider the complex signal. 


z{t) = aocos{uot) + jaaosm{uot) 


(62) 


and the spectrnm given by (60) which is not one-sided. However in terms of a single AM-FM 
component, we can have 


z{t) = 


(63) 


where 


ao{t) = \Jal cos2(a;ot) 


-f- sin" 


(wo«) 


(64) 


and 


6*0 (t) = arctan 


aoo sin(a;ot) 
Oo cos(a;ot) 


(65) 


The IF, ujoit) 


d 

dt 


6 * 0 (t) is strictly positive for any positive choice of a. Thns the associated 


Hilbert spectrnm is one-sided even thongh the spectrnm is two-sided. In other words, saying 
we assnme positive IF is not the same as saying we assnme a one-sided spectrnm. Althongh 
signals can be designed snch that Z{u), for practical pnrposes, has a one-sided spectrnm, 
e.g. commnnications signals, this does not mean that natnrally-occnrring signals have, in 
general, a one-sided spectrnm. Therefore, making the one-sided spectrnm assnmption (and 
assnming HC) may lead to incorrect model parameters. 

Another incorrect assnmption that is often made is that if the signal is narrowband, the 
HT will yield a meaningfnl complex extension. Unfortnnately, this is not the case becanse 
narrowband signals exist that do not have Hermitian symmetry and hence nse of the HT will 
yield incorrect resnlts. This is demonstrated by (60) which conld be considered narrowband, 
yet use of the HT to extend its real observation would not yield the correct latent signal due 
to the absence of HC. 

By changing from SHCs to AM-FM components and relaxing the assumption of HC, we 
have generalized the spectrum to a Hilbert spectrum and effectively changed the dehnition 
of IF to be component-dependent. By reverting back to Carson’s dehnition and not using 
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SHCs in the analysis, we move away from Gabor and Ville’s specialized definition back 
towards the generalized definition of IF. 

While using AM-FM components to define IF, which can change in time, we have to 
consider the possibility that a component’s IF changes sign. Such a sign change is not possi¬ 
ble under the assumption of SHCs which have constant IF. In this work, we have arbitrarily 
assumed that all components must have non-negative IF for all t, although there may be 
some signal classes for which this is not true. In these classes, the AM-FM parameters may 
not match the underlying signal model parameters, although the superposition will yield the 
correct real signal. This is illustrated in Fig. 7 where a component of z{t) with associated IF 
u{t) and the component of z*{t) with associated IF —oj{t) cannot be separated by assum¬ 
ing non-negative IF because each has both positive and negative instantaneous frequencies. 
Therefore in these cases, we need to relax the the assumption of non-negative IF at some 
time instances in order to properly estimate the latent signal. 


6. Examples using the AM—FM Model 

The example in Part I illustrated solutions to the LSA problem where we arrived at a 
single pair of lA/IF parameters, p{t) and G(f), for the latent signal. In this section, we 
illustrate HSA consisting of a set of lA/IF pairs, {ak(t), cokit)} for the signal, one pair for 
each latent component. This is illustrated by the set diagram in Fig. 3. In our closed-form 
solutions, we choose assumptions which lead to various decompositions and interpretations. 
For convenience, it is assumed that the phase reference <fk = 0, where possible. 


6.1. Periodic Triangle Waveform Example 

As an example of HSA, we consider the triangle waveform from Part I (24) and the 
sinusoidal FM signal. 


6.1.1. Analysis assuming simple harmonic components 

If we assume the component in (36) has the form given in (45), i.e. a SHC, then the 
AM-FM model can be obtained using Fourier analysis and the triangle waveform may be 
expressed with an inhnite number of components 




^j(2k+l)uJot 


] 


where the lA is given by 


the IF is given by 


ak{t) = 


8A 1 

TT^ {2k 1)^ ’ 


uJk{t) = {2k + l)a;o, 

and uq is a constant—conventionally known as the fundamental frequency. 


( 66 ) 


(67) 

( 68 ) 


25 




6.1.2. Analysis assuming a single AM-FM component with harmonic correspondence 

The solution in (66) may be collapsed into a single, wideband AM-FM component (K = 
1) by rewriting as in Part I, in terms of the AM-FM model as 

x{t) = 3f? (69) 

where the lA ao(t) is given by p(t) in (26) and the IF oj(t) is given by Cl(t) in (27). 

6.1.3. Analysis assuming a single AM component 

If we assume a single component {K = 1) with constant frequency u!o{t) = uq, the 
triangle waveform may be expressed as in Part I, in terms of the AM-FM model as 

x(t) = 3fJ{ao(t)e^‘"“'} (70) 

where the lA ao{t) is given by p{t) in (30). 

6.1.4. Analysis assuming a single FM component 

If we assume a single component (K = 1) with constant amplitude ao{t) = A, the triangle 
waveform may be expressed as in Part I, in terms of the AM-FM model as 

a:(t) = (71) 

where the IF oJo(t) is given by fl(t) in (32). 


6.2. Sinusoidal FM Example 

A sinusoidal FM signal with carrier frequency cuc, modulation frequency Um, and constant 
B, is expressed as 

x{t) = | _ (^72) 


6.2.1. Analysis assuming a single FM component 

We recognize that (72) is already in the form of the AM-FM model, if we assume a single 


component (K = 1) with lA 

ao(t) = 1. 

(73) 

This leads to the IF ^ 

Uo{t) =Uc+ —Bsm{Ujnt). 

dt 

(74) 


6.2.2. Analysis assuming simple harmonic components 

Alternatively, the sinusoidal FM signal can be expressed in terms of the AM-FM model 


as 


a;(t) = 3?< ^ .Jk{2nB/ujm)e 




(75) 


, k=—oo 
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where Jk{ ) denotes the fcth-order Bessel function of the hrst kind [19]. This yields an inhnite 
number of components with lA given by 

ak{t) = Jk{‘27iB/um) (76) 

and IF given by 

^kit) ku^rfYi. (77) 

This is an example of a signal class where the assumption of non-negative IF components 
in (41) must be relaxed, as a consequence of the summation on k becoming double-sided. 
For most practical choices of B and cjm, the Jk{27rB/um) associated with the negative 
frequencies will be vanishingly small, therefore the assumption of non-negative IF, is for 
practical purposes, true. 

6.3. Remarks 

As demonstrated above, we obtain different parameterizations of the same signal by 
changing assumptions made during analysis. Some parameterizations, such as those obtained 
using Fourier analysis and the HT, are closely related due to similar assumptions. However, 
there exist many parameterizations each with different interpretations and without other 
information, we cannot say that one particular parameterization is any better or worse 
than any other. This highlights the importance of considering the hidden assumptions when 
using any particular method (FT, HT, STFT, HHT, etc.) and interpreting meaning from the 
parameters. If assumptions in analysis do not match the underlying signal model, erroneous 
interpretations may be made. 

The previous solutions make use of different quadratures each corresponding to different 
signal models. This demonstrates the non-uniqueness of the quadrature for the AM-FM 
component, and hence our view of the quadrature as a free parameter. Also, note that the 
solutions in (66), (69), (70), etc., each utilize a complex extension that is implied by the 
particular assumptions made, rather than utilizing an extension based on a single rigorously- 
dehned procedure, like the HT. In fact, for any practically-chosen real signals x{t) and y{t), 
there exist assumptions leading to instantaneous parameters in which y{t) can be viewed as 
the quadrature of x{t). 

The AM-FM model leads to exact solutions for a{t) and uj{t), so it might appear to 
violate the uncertainty principle, i.e. exceed the lower limit of the time-bandwidth prod¬ 
uct. However, when AM-FM modeling is viewed as a quantum mechanics problem, our 
casting of the problem is fundamentally different than Gabor’s. A comparison of the formal 
correspondences between quantum mechanics and short-time Fourier analysis and quantum 
mechanics and Hilbert spectral analysis is given in Table 2. In Gabor’s casting, time and 
frequency are the complementary variables, i.e. (position, momentum)<H-(f, u) [1, 88] while 
for the AM-FM model, the real and quadrature signals are the complementary variables, 
i.e. (position, momentum)<H-(a;(f), y{t)). Therefore, all uncertainty arises because only the 
real signal is observed. We believe that our casting is not only more appropriate, but also 
more useful and powerful. 
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Table 2: The formal correspondences between quantum mechanics and short-time Fourier analysis (p. 197 
in [1]) (columns 2 and 1) and quantum mechanics and Hilbert spectral analysis (columns 2 and 3). 


Fourier Analysis 


Quantum Mechanics 


Hilbert Spectral Analysis 

Description 

Symbol 


Description 

Symbol 


Description 

Symbol 

Time 

t 

■fA 

Position 

<? 

■fA 

Real Signal 

x{t) 

Frequency 

U) 

•H- 

Momentum 

p 

•H- 

Imaginary Signal 

y{t) 

no correspondence 


•H- 

Time 

t 

•H- 

Time 

t 

Signal 

x{t) 

■fA 

Wave Function 


■fA 

Latent Signal 

z{t) 

Uncertainty 

ST> 1 

■fA 

Uncertainty 

O'pCTq ^ 2^ 

■fA 

Uncertainty 

Instantaneous Frequency 

x{t) = 5R{2(t)} 

= it argfbfch)} 


7. Proposed Visualization of the Hilbert Spectrum 

The ability to visualize and interpret model parameters is key to the adoption of any 
analysis method. Often complex AM-FM signals are plotted as a series of real 2-D plots, 
i.e. Sk{t) vs. t, which for AM-FM components, would provide little insight into the underlying 
signal model. Alternatively, Argand diagrams may be used for AM-FM signal visualization, 
however, drawbacks include the quadrature signal possibly having no assigned meaning and 
the revolution rate not intuitively displayed. 

A more appropriate visualization plots the Hilbert spectrum in its entirety. Unfortu¬ 
nately, plots of the Hilbert spectrum are often crudely discretized because a clear distinction 
between instantaneous parameters and spectral parameters is not made [42]. We propose a 
method for visualizing the Hilbert spectrum, which is both complete, intuitive, and avoids 
the coarse discretization. The proposed visualization can be considered a (pseudo-) phase- 
space plot because every degree of freedom or parameter of each AM-FM component is 
represented. 

By plotting uJk{t) vs. Sk{t) vs. t as a line in a 3-D space and coloring the line with respect 
to \ak{t) \ for each component, the simultaneous visualization of multiple parameters for each 
component is possible. Further, orthographic projections yield common plots: the time-real 
plane (the real signal waveforms), the time-frequency plane (Hilbert spectrum), and the 
real-frequency plane (analogous to the Fourier magnitude spectrum). We have found it 
benehcial to interpret each component as an “illuminated” particle moving in 3-D space. 
Each particle’s motion in time and frequency is governed by The proposed method 

allows one to visualize the assumed underlying signal model. 

To illustrate visualization of the Hilbert spectrum, we plot the examples in Section 6. 
However, the sophisticated nature of the proposed 3-D visualization of the Hilbert spectrum 
is not well-accommodated with paper media, and as a result we provide associated matlab 
functions and additional visualizations online [89]. It is our preference that Hilbert spectra 
not be visualized with paper media. In the plots, we have utilized a perceptually-motivated 
colormap in order to improve interpretation over other colormaps [90, 91]. 

For the triangle waveform example. Figs. 8(a),(c),(e) illustrate the Hilbert spectrum 
plots for the assumptions of SHC, single AM-FM component with HC, and single AM 
component, respectively. Fig. 8(a) shows the hrst three SHCs [constant amplitude and 
constant frequency in (67)] at integer multiples of a fundamental frequency, hOvr rads/s. 

28 








Fig. 8(c) shows the single AM-FM component with HC, where we see line color variation 
indicating a time-varying lA in (26) and a clear time-varying IF in (27). Fig. 8(e) shows 
the single AM component, where we see a constant IF and color variation indicating the 
time-varying lA in (30). Figs. 8(b),(d),(f) show the corresponding time-frequency planes of 
Figs. 8(a),(c),(e) by projecting out the Sk(t) dimension. 

For the sinusoidal FM example, Figs. 9(a) and (c) illustrate the Hilbert spectrum plots 
for the assumptions of a single FM component and SHCs respectively. Fig. 9(a) shows the 
constant lA in (73) indicated by a constant line color and time varying IF in (74). Fig. 9(c) 
shows SHCs at integer multiples of a fundamental frequency 47r rads/s where each component 
has constant lA in (76). Figs. 9(b) and (d) show the corresponding time-frequency planes 
of Figs. 9(a) and (c) by projecting out the Sfc(t) dimension. 

8. Summary 

In Part H of this paper, we have presented theory for the Hilbert spectrum as a gener¬ 
alized LSA problem. In the general problem, we seek a representation of z{t) consisting of 
a superposition of latent signals, i.e. multicomponent model consisting of a superposition 
of complex AM-FM components. We have used the AM-FM model to dehne the Hilbert 
spectrum as parameterized by a set of lA/IF pairs (and phase references) each associated 
with the components. In the LSA problem, relaxation of the HC condition allowed freedom 
in the choice of the signal’s quadrature and admitted many solutions for the latent signal. 
In the HSA problem, relaxation of the HC condition allows freedom in the choice of each 
component’s quadrature thereby allowing even greater freedom in the construction of the 
signal’s model. We illustrated how assumptions on the form of the component, i.e. constant 
amplitude and constant frequency, time-varying amplitude and constant frequency, constant 
amplitude and time-varying frequency, etc., lead to familiar specializations of the model. 

We discussed the frequency domain view of LSA including the implications of relaxing 
HC. From there we discussed the Hilbert spectrum view of LSA, specihcally the implications 
of using a general AM-FM component with relaxed HC and assumed positive IF. Closed form 
HSA examples were provided in which we assume various forms of the AM-FM component 
and determined the unique corresponding instantaneous parameterization in terms of lA 
and IF. A novel 3D visualization of the Hilbert spectrum was proposed by plotting u{t) 
vs. s{t) vs. t and coloring with respect to |a(t)|. 

Finally, we discussed the analogy of the Hilbert spectrum to quantum mechanics in which 
our casting of time-frequency analysis is fundamentally different than Gabor’s casting, where 
the uncertainty in our framework is in the quadrature signal and not the frequency variable. 
This provides a new and powerful framework for nonstationary signal analysis. 
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Figure 6: Illustrations of when Gabor’s QM (a) can distinguish Z{uj) = Zi{uj) + Z2{(jj) from Z*(—uj) = 
Z^(—uj) + Z^iy—oj) because Z{ijj) has HG, (b) cannot distinguish Z{ijj) = Zi^oj) + ^ 2 ( 0 ;) from Z*{—uj) = 
ZJ‘(—w) + Z^i—uj) because the latent spectrum is two-sided and Hermitian symmetry is imposed by the 
real observation, and (c) is incorrect because the structure of the latent spectrum along with the Hermitian 
symmetry imposed by the real observation has concealed terms. 
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Figure 7: Illustration of when the assumption of non-negative IF cannot distinguish z{t) with associated IF 
oj{t) (—) from z*{t) with associated IF —oj{t) (—) because each has both positive and negative values of IF 
at some time instances. 
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Figure 8: Hilbert spectrum for the triangular waveform x{t) with wg = SOtt rads/s for the assumptions 
of (a) SHCs, (c) single AM-FM component with HC, and (e) single AM component. The corresponding 
time-frequency planes, obtained by projecting out the Sk(t) dimension, are shown in (b),(d),(f). 
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(c) (d) 


Figure 9: Hilbert spectrum for the sinusoidal FM waveform x{t) with uJc = IIOtt rads/s, uJm = 47r rads/s, and 
B = 2b for the assumptions of (a) a single FM component and (c) SHCs. The corresponding time-frequency 
planes, obtained by projecting out the Sk{t) dimension, are shown in (b),(d). 
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Part III: Numerical HSA Assuming Intrinsic Mode Functions 

[Tjhere are some crucial restrictions of the Fourier spectral analysis; the system 
must be linear; and the data must be strictly periodic or stationary; otherwise, 
the resulting spectrum will make little physical sense.—Huang [42] 

In Part II, the AM-FM signal model and HSA theory was presented, analysis was per¬ 
formed on example signals, and a visualization of the Hilbert spectrum was proposed. In 
Part HI, numerical algorithms are given for estimating the lAs and IFs of the components in 
the AM-FM model where now the component is assumed to be an Intrinsic Mode Function 
(IMF). These algorithms hrst decompose the signal into IMFs using an improved version of 
Huang’s original Empirical Mode Decomposition (EMD) algorithm and second, demodulate 
the IMFs to obtain the instantaneous parameters. Unlike previous studies, close attention 
is paid to the assumptions made in the dehnition of the IMF which are carried forward 
to the demodulation step, thereby avoiding any ambiguity associated with obtaining the 
instantaneous parameters. We begin with a comprehensive review of EMD and several vari¬ 
ations, and propose an algorithm for the computation of the Hilbert spectrum assuming 
IMF components. It is important to note that while IMFs can be considered latent AM-FM 
components, there are other classes of AM-FM components that are not IMFs as illustrated 
in Fig. 5. Examples using the proposed algorithm are provided that highlight alternative 
decompositions compared to traditional Fourier analysis and demonstrates the advantages 
of using the HSA framework. 

1. Introduction 

In Part H of this paper, we dehned the Hilbert spectrum using the AM-FM model in 
(35) where the AM-FM component is dehned in (36) and parameterized by the lA ak{t), IF 
cokit), and phase reference (pk- We assume a real observation, x(t) of a latent signal, z(t). 

In Part HI of this paper, we turn our attention to numerical computation of instantaneous 
parameters of the AM-FM model. Many methods for computing the parameters of an AM- 
FM model already exist, each corresponding to a specihc set of assumptions. As discussed 
in Part H, a FT corresponds to the assumption of SHCs and a STFT corresponds to the 
assumption of AM components. Other AM-FM models have been proposed with alternative 
assumptions, however, common to these are restrictions which limit the utility of the model. 

Practical estimation of the instantaneous parameters of the AM-FM model in (35) is 
a two-step process. First, the signal must be decomposed into a set of latent AM-FM 
components and second, the instantaneous parameters {ak{t),Uk{t)} of each component 
must be estimated. To date, the most hexible decomposition method for AM-FM modeling 
is Huang’s EMD and its variations [42] . In contrast to most time-frequency analysis methods, 
EMD makes no assumption of constant amplitude or frequency on the component and has 
mild bandlimiting assumptions. In [42], Huang proposed the original EMD algorithm which 
sequentially determines a set of AM-FM components, i.e. IMFs via an iterative sifting 
algorithm. The Ensemble Empirical Mode Decomposition (EEMD) [92] and tone masking 
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[93] introduced ensemble averaging in order to address the mode mixing problem. The 
complete EEMD was proposed to address some of the undesirable features of EEMD by 
averaging at the component-level as each component is estimated rather than averaging at 
the conclusion of EMD [94] . Several improvements to the sifting algorithm have also been 
proposed including those by Rato [95]. These improvements to the original EMD will be 
more fully described in this part. 

In the second step, instantaneous parameters must be estimated through demodulation of 
the IMFs. Despite the numerous attempts to demodulate IMFs, most of the proposed meth¬ 
ods have fallen short because the assumptions made during decomposition have not been 
maintained during demodulation. In this paper, we identify and pay close attention to these 
assumptions as we develop a demodulation technique which adheres to the assumptions. We 
point out that Rato proposed an AM-FM demodulation procedure in which the lA estima¬ 
tion was consistent with the assumptions but the IF estimation was not [95]. Also, Huang 
has examined numerous demodulation methods, including an iterative normalization to ob¬ 
tain an FM signal which is then demodulated to estimate IF using an arctan approach that 
unfortunately suffers from numerical instability [96] . The iterative normalization method is 
consistent with the decomposition assumptions. Utilizing both Rato’s AM estimation and 
Huang’s iterative normalization procedure, we propose a mathematically equivalent method 
to get IF from the FM signal that is more numerically stable. We then incorporate the 
proposed demodulation and EMD into a single HSA-IMF algorithm which gives very good 
estimates for the lA and IF parameters of the AM-FM model. 

Part HI of this paper is organized as follows. In Section 2, we review the EMD algorithm 
and its variations and improvements by other researchers. By carefully noting the assump¬ 
tions made in EMD, we propose demodulation of the resulting IMFs and present a complete 
HSA-IMF algorithm composed of EMD and demodulation. In Section 3, we provide two 
sets of examples using the HSA-IMF algorithm. The first set of examples uses synthetic 
signals which illustrate how EMD can estimate the parameters of the underlying, assumed 
AM-FM signal model. The second set of examples uses real-world audio signals including 
speech which further illustrates AM-FM model parameterizations with alternate interpreta¬ 
tions. These examples clearly demonstrate the alternate decomposition and interpretation 
offered by HSA. In Section 4, we discuss other issues related to HSA and in Section 5, we 
discuss future research. In Section 6 we summarize Part HI. 

2. Empirical Mode Decomposition 

EMD consists of an iterative procedure for decomposing a signal into a set of IMFs, 
{<Pk(t)} [42]. In [42], the definition of an IMF is any signal that satisfies two conditions: 

Cl: In the whole signal segment, the number of extrema and the number of zero crossings 
must be either equal or differ at most by one. 

C2: At any point the mean value of the envelope, defined by the local maxima and the 
envelope defined by the local minima, is zero. 
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In the context of LSA and AM-FM modeling, an IMF can be thought of as a latent compo¬ 
nent, 

ipkit) = (78) 

Thus we view an IMF as an inherently comp lex-valued component, which is not the conven¬ 
tional interpretation. Furthermore, we argue that the definition of an IMF forces a unique 
quadrature that in general does not equal ^{{ipkit)} as will be discussed in Section 2.4. 

Unlike traditional methods in time-frequency analysis, EMD is dehned by an algorithm 
rather than transform theory, which has both advantages and disadvantages. One advantage 
is that EMD utilizes a component that is far less restrictive than other time-frequency meth¬ 
ods [1]. One disadvantage is that EMD is understood primarily through experimentation 
[93]. Empirical experiments using white noise have shown EMD to act as a dyadic filter 
bank [97, 98, 99, 92, 100]. Using fractional Gaussian noise as a model for broadband noise, 
it has been shown that the built-in adaptivity of EMD makes it behave spontaneously as a 
‘wavelet-like’ hlter, i.e. the result of sifting is a ‘detail’ and a ‘trend’ [98]. 

Efforts have also been made to replace the sifting algorithm with alternate formula¬ 
tions which are more mathematically grounded, such as techniques based on optimization 
[101, 102, 103], machine learning [104], PDEs [105, 106, 107, 108, 109, 110], and Fourier anal¬ 
ysis [111]. We also point out other research on the EMD algorithm including the multivariate 
EMD [112, 113, 114, 115, 100, 116, 117] and the complex EMD [118]. As a final note, not 
all variations of EMD utilize the IMF component or the proposed AM-FM model and thus 
cannot be considered as a form of HSA. Examples include the Hilbert variational decom¬ 
position [119, 81, 66], the time-dependent intrinsic correlation [120], and synchrosqueezed 
wavelet transforms [121, 122, 41]. 

2.1. The Original EMD Algorithm 

The original EMD and sifting algorithms proposed by Huang [42] are listed in Algorithms 
1 and 2. The purpose of the sifting algorithm is to iteratively identify and remove the trend 
from the signal, effectively acting as a high pass hlter. Step 5 of Algorithm 1 removes the 
high frequency component (fkit), estimated during the sifting process. The process is then 
repeated to remove additional IMFs from the signal if they exist. The resulting decomposi¬ 
tion is complete and sparse [42, 123, 103]. EMD is formulated with continuous-time signals, 
however, in practice, EMD is applied to discrete-time signals which may result in errant 
decompositions. The effects of sampling in the context of EMD have been considered by 
Rilling and it is generally recommended to oversample but not resample before application 
of EMD, so that EMD effectively behaves like a continuous operator [124]. 
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Algorithm 1 Empirical Mode Decomposition 
1; procedure {^Pk{t)} = EMD( x{t) ) 

2; initialize: k = 0 and x_i(t) = x{t) 

3: while Xk-\ (t) 4^ 0 and Xk-\ it) is not monotonic do 

4: (Pfc(t) = SIFT( Xk-i{t) ) 

5: Xk{t) = Xk-i{t) - ^pk{t) 

6: k k + 1 

7: end while 

8: (pfe(t) = Xfc-i(t) 

9: end procedure 


Algorithm 2 Sifting Algorithm 
1: procedure (p(t) = SIFT( r{t) ) 

2: initialize: e{t) 4 0 

3: while e{t) ^ 0 do 

4: find all local maxima: Up = r{tp), p = 1,2,... 

5: find all local minima: Ig = r{tq), q = 1,2,... 

6: interpolate: u{t) = CublicSpline({fp, Mp}) 

7: interpolate: l{t) = CublicSpline({fg,/g}) 

8: e{t) = [u{t)+ l{t)]/2. 

9: r(t) r(t) — e(t). 

10: end while 

11: ip{t)=r{t) 

12: end procedure 


2.2. Improving the Sifting Algorithm 

If EMD is viewed as an AM-FM decomposition technique, then the sifting algorithm 
is an iterative way of removing the asymmetry between the upper and lower envelopes in 
order to transform the input r{t) into an IMF [95]. By doing so, low frequency content is 
discarded at every sifting iteration, effectively making the sifting algorithm behave as a high 
frequency hlter or high frequency component tracker. Due to the doubly iterative nature of 
EMD and termination conditions, numerical imprecision and differing implementations can 
lead to very different IMFs. To achieve consistency when using EMD, Rato proposes the 
following constraints [95]: 

Scale: IMFs should scale with the signal. 

Bias: Any signal bias, should only be reflected in the trend. 

Identity: The EMD of an IMF should be the IMF itself.® 


®As Rato points out, this constraint may be relaxed in some cases [95]. 
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Time-reversal: Time-reversal of the signal should time-reverse the IMFs. 

Several improvements have been made to the sifting algorithm (Algorithm 2) to improve 
the decomposition accuracy [125, 95]: 

1. Improving stop criterion robustness in Step 3 

2. Providing a consistent method for identihcation of extrema in Steps 4 and 5 

3. Addressing interpolation end effects in Steps 6 and 7 

4. Scaling the mean envelope removal in Step 9 

Although several stopping criteria have been proposed [126, 127, 95], Rato suggested the use 
of a resolution factor which is the ratio between the energy of the signal at the beginning 
of sifting r{t) and the energy of the average of the envelopes e{t). If this ratio increases 
above a predetermined threshold, then the IMF computation terminates. We have found 
that Rato’s use of a parabolic interpolator [95] to identify the extrema in Steps 4 and 5 is 
convenient because it uses as few as three samples and it avoids the classihcatory function 
to determine if a particular sample is a maxima, minima, or neither. End effects appear 
due to the fact that a given interpolator may not be a good extrapolator [95]. In order to 
deal with this problem, Rato suggested to insert artihcial minima and maxima in order to 
control the behavior of the interpolator. In addition, the mean envelope removal is scaled 
in Step 9, 

r{t) r{t) — ae{t) (79) 

where the step-size 0 < a < 1 increases the number of sifting iterations but improves 
stability and robustness of the resulting IMFs [95]. 

We illustrate a single iteration of the sifting algorithm in Fig. 10. In Fig. 10(a), we 
show two unknown components, (foit) ( — ) and (fiit) (— ) and in Fig. 10(b), we show the 
signal under analysis r(t) = (po(t) -|- (pi(t) (—) which is the input to the sifting algorithm. 
Figs. 10(c) and (d) show the location and interpolation of the extrema as given in Steps 
4-7. Interpolations lead to estimates of the upper envelope u(t) (-^) and lower envelope 
l(t) (-^) of r(t). The average of the upper and lower envelopes e(t) ~ </?i(t) is shown in 
Fig. 10(e). At the end of the hrst iteration of sifting, r{t) — e{t) ~ ‘^o(^) and is shown in 
Fig. 10(f). 

2.3. Improving the EMD Algorithm 

The major problem in the EMD algorithm is mode mixing, which is dehned as a sin¬ 
gle IMF either consisting of components of disparate scales or components of similar scale 
residing in the same IMFs [92]. Mode mixing is a consequence of signal intermittency, or 
more specihcally relative component intermittency. As a result, the particular component(s) 
tracked by the sifting algorithm in a particular IMF at any instant may change as intermit¬ 
tent components begin or end [93]. This is illustrated in Fig. 11(a), where in the left part 
of the hgure, we show components of disparate scales being in the same IMF denoted by O, 
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while in the center part of the hgnre we show two components of similar scale in the same 
IMF denoted by □. The ability of EMD to resolve two components considering both the 
relative lAs and IFs of components, was examined and qnantihed by Rilling [128]. However 
as was noted, resolving closely-spaced components may not be the nltimate goal, provided 
that the decomposition is snitably matched to some meaningfnl interpretation [128]. 

Two commonly nsed methods of mitigating mode mixing are EEMD [92] and tone mask¬ 
ing [93]. EEMD (given in Algorithm 3, where E[-] denotes the expectation) ntilizes zero-mean 
white noise, to pertnrb the signal so a component may be tracked properly over an 

ensemble average. As the illnstration in Fig. 11(b) shows, noise can be nsed to assist the 
sifting algorithm. Inserting noise with high enongh power gives the sifting algorithm some¬ 
thing to track when the highest freqnency component is intermittent, then vanishes in the 
ensemble average. Althongh injecting noise can help to track components properly, a care- 
fnlly designed masking signal can resnlt in better performance. A sitnation where a carefnlly 
designed masking signal may be benehcial is illustrated in Fig. 11(a) denoted by □. 


Algorithm 3 Ensemble Empirical Mode Decomposition 
1: procedure {(pfc(^)} = EEMD( x{t) ) 

2: initialize: I is the number of trials 

3: ^^\t) = EMD ( xit) + n;«(t) ) , * = 1,..., J 

4: =E[ip^^\t)] 

5: end procedure 


EEMD is not without its disadvantages as it is more computationally complex, loses 
the perfect reconstruction property, propagates IMF estimation error, results in inconsistent 
numbers of IMFs across the trials, and the resulting set of averaged IMFs {<fk{t)} are not 
necessarily IMFs [94]. Torres proposed the “complete EEMD” given in Algorithm 4 to 
address some of these issues [94]. Complete EEMD dehnes a procedure EMDfc(-) which 
returns the /cth IMF obtained using EMD [94]. This method of EEMD requires fewer sifting 
iterations, a smaller ensemble size, and recovers the completeness property of the original 
EMD algorithm to within the numerical precision of the computer [94]. 


39 






Algorithm 4 Complete EEMD 
1; procedure = CEEMD( x{t) ) 

2; initialize: A; = 1, is a SNR factor, 

and I is the number of trials 
1 ^ 

3: ipoit) = - 5^SIFT( xit) + /3o^»(f) ) 

i=l 

4: Xoit) = x(t) - (po(t) 

5: while Xk-i(t) ^ 0 and Xk-iit) is not monotonic do 

1 ^ 

6: (pfc(t) = - 5^SIFT( Xk-i{t) + /3fcEMDfc( ) 

i=l 

7: Xkit) = Xk-iit) - (pkit) 

8: k k + 1 

9: end while 

10 : (pk{t) = Xk-i{t) 

11: end procedure 


Rather than using a noise signal, a deterministic signal v{t) can also be used as a per¬ 
turbation and then removed after IMF estimation. The “tone masking” technique is given 
in Algorithm 5 and can be used, for example, in place of Step 4 in Algorithm 1 [93]. Ad¬ 
vanced forms of tone masking have been proposed and are termed Signal-Assisted EMD 
(SA-EMD) [99, 129, 130, 131]. These methods have additional advantages in their ability 
to track closely-spaced components, but require careful selection of the masking signal. 


Algorithm 5 Tone Masking 
1: procedure (p[t) = TM( x{t),v{t) ) 

2: = x{t) + v{t) and x^~\t) = x{t) — v{t) 

3: = SIFT( ) and = SIFT( x^-\t) ) 

4: m 

5: end procedure 


2 .J 1 .. IMF Demodulation 

In order to obtain the lA/IF parameters in the AM-FM model, we are required to 
demodulate the IMFs returned by EMD. One approach, used by Huang, is to apply the HT 
to each IMF in order to obtain estimates of lA and IF. This analysis is often referred to as the 
Hilbert-Huang Transform (HHT) [42]. Using the HT for IMF demodulation is inappropriate 
because as argued, the HT assumes SHCs and HC—which is not true since an IMF is in 
general an AM-FM component without HC. A similar argument can be made of many of 
the other demodulation methods that have been considered [96]. A second approach, as 
advocated in Part H, is to let the assumed form of the AM-FM component imply a complex 
extension. When viewed in the context of LSA, the definition of the IMF given in Section 
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2 forces a unique complex extension to the IMF—justifying our view of the IMF as a latent 
component specihed by a real signal. 

Inherent in C2 of the IMF dehnition, are the following assumptions: 

Al: a{tp) = |s(tp)|, where \s{tp)\ are the extrema of s{t). 

A2: a{t),t ^ {tp} is inferred by cubic spline interpolation. 

The hrst assumption, Al can be viewed as quadrature forced to zero at {tp}. To see this, 
note that from (36), we have 

|a(t)| = (80) 

and thus, <j(tp) = 0 and a(tp) = |s(tp)|. Al also implies non-negativity of a(tp), however, 
non-negativity of a{t) is not guaranteed for all t due to the interpolation. 

The second assumption, A2 can be viewed as a relative bandlimiting condition on a{t) 
controlled by two factors: the choice of interpolator and the density of extrema. To see 
this, note that a{t) between extrema of s{t) is dehned by an interpolator. This implies that 
a{t) will only have as much variation between extrema as the interpolator allows. However, 
with dense extrema (high IF) much less restrictive constraints have been imposed on a{t) 
than if extrema are sparse (low IF)—thus our view as a relative bandlimiting condition on 
a{t). Also note, using an interpolator other than the cubic spline is likely to change the lA, 
effectively changing the resulting IMFs [95, 132]. 

In the IMF dehnition. Cl requires that the number of extrema and the number of zero 
crossings must be either equal or differ at most by one. A general AM-FM component may 
not satisfy this condition, for example as a result of sign reversal of u:{t) or in cases where 
Al and A2 are not satished. Assumption of positive IF on the AM-FM component in Part 
II, eliminates the possibility of sign reversal of the IF but may not be true for all signal 
classes. 

Rato proposed an AM demodulation approach given in Algorithm 6, that is consistent 
with the decomposition assumptions in Algorithm 2 [95]. Starting with an IMF estimate 
we obtain an estimate for lA a{t) which can then be used to estimate the IF via the 
FM signal 

■Sfm(^) = (81) 

a(t) 

However, the estimate in (81) may result in |sfm(^)| > 1- Thus, Huang proposed an iterative 
normalization procedure, given in Algorithm 7, which removes the AM from the signal to 
obtain a more accurate sfm(^) [96]. Although the iterative procedure improves demodula¬ 
tion accuracy, we noticed that it can be susceptible to oscillating artifacts introduced by 
overhtting of the cubic spline interpolator. As a result, we have found that these artifacts 
can be minimized by limiting the number of iterations to three. 
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Algorithm 6 IMF lA Estimation 
1; procedure d{t) = IAest( (p{t) ) 

2: r{t) = \(p{t)\ 

3: find all local maxima: Up = r{tp), p = 1,2,... 

4: interpolate: d(t) = CublicSpline({tp, Mp}) 

5: end procedure 


Algorithm 7 Obtaining a real FM signal from an IMF 
1: procedure sfm(^) = iterAMremoval( (p{t) ) 

2: initialize: g{t) = (p{t), b{t) ^ 1, and n = 0 

3: while h{t) ^ 1 and n < 3 do 

4: hit) = IAest( git) ) 

5 : g{t) ^ 9{t)/K'i) 

6: n •(— n + 1 

7: end while 

8: SfmW = 9{t) 

9: end procedure 


We can directly obtain two estimates of the IF, ±a)(t) by substituting ait) = 1 in 
(36b), sit) = sfm(^) in (36c), and then computing (41). Direct FM demodulation is not 
straightforward because different approaches exist for obtaining the IF from spuit), that 
although are mathematically equivalent, may differ in numerical stability [96]. 

We propose a FM demodulation method to address the numerical stability issues. We 
begin by estimating the quadrature in (36c) as 


dFM(t) = -sgn 


d 

dt 


Spuit) 




(82) 


where —sgn [^SFM(t)] is required to obtain an appropriate four quadrant estimate with 
assumed positive IF. In Fig. 12 , we have provided an illustration (animated in the web ver¬ 
sion) to assist the reader with understanding of how we arrive at (82). The computationally 
unstable points, {to} occur where (Tfm(^o) = 0, thus we can replace a small range around 
these points (to — e, to + e) with interpolated values. Then 

6{t) = arg [sFM(t) + jo-puit)] (83) 


and the IF is obtained with (41). We can optionally smooth the resulting IF estimate. Our 
complete IMF demodulation algorithm is listed in Algorithm 8. 
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Figure 12: In this figure (animated in web version), SFM(i), d'FM(i) are represented by the blue and green 
vectors, respectively. The amplitude-normalized 'ipFAiit), is represented by the red vector. The magnitude 
of apM{t) is easily calculated. The sign of (TFM(i) is obtained as follows. We note that when SFM(i) is 
decreasing (blue vector moves left), i3'fm( 0 is always positive (green vector is in the upper half plane) and 
when spuit) is increasing (blue vector moves right), (TFM(i) is always decreasing. Thus by reversing the sign 
of the derivative of SFM(i) we obtain the sign of ctfm(0- 


We note to the reader that the proposed IMF demodulation procedure, which does not 
involve the HT or FT, can in some cases return the same results but in general this is 
not true. For example, demodulation of cos(a;ot) with either the IMF or HT demodulation 
returns a SHC with HC. On the other hand, HT demodulation of the triangle waveform 
returns a single AM-FM component with HC as given in (27), while IMF demodulation 
returns IF given by (32) and constant lA. 

2.5. Proposed Algorithm for HSA Assuming IMFs 

The various recommendations and improvements listed this section appear to have never 
been combined and investigated together. Therefore, we propose to incorporate the most 
desirable features of the complete EEMD and tone masking for addressing the mode-mixing 
problem, Rato’s improvements to the sifting algorithm, and our proposed demodulation 
method, into a single algorithm given in Algorithm 9. 
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Algorithm 8 IMF demodulation 


procedure [a{t),u{t)] = IMFdemod( (p{t) ) 
a{t) = IAest( (p{t) ) 

’Sfm(^) = iterAMremoval( (p{t) ) 
d ^ / \ 


4: apuit) = -sgn 


dt 


12 


’FM 


(t) 


5: Find {to} such that (TFM(to) = 0 

6; For each to, replace (dFM(^o — e), dFM(^o + e)) 
with interpolation 

7: a)(t) = —arg [sFM(t) + jdFM(t)] 

8 : end procedure 


Algorithm 9 HSA-IMF Algorithm 
1 : procedure {(pfc(t), afc(t), a)fc(t)} = HSA-IMF( x{t) ) 

2: initialize: a;_i(t) = x(t), k = 0, (3k is a. SNR factor, e: is an energy threshold, and / 

is the number of trials 
3: while J \xk-i{t)\'^dt > e 

and Xk-i{t) is not monotonic do 
1 ^ 

4 : Mt) = ) 

5: [afc(t),a;fc(t)] = IMFdemod( ipk{t) ) 

6 : Xk(t) = Xk-i(t) - (fkit) 

7: k k + I 

8 : end while 

9: (pfe(t) = Xfc-i(t) 

10: end procedure 


Although the masking signal is usually a carefully chosen real-valued AM-FM 

signal [129, 130], we choose as a white, Gaussian, lowpass-hltered noise with cutoff 

frequency below the amplitude-weighted IF [133] of the previously found IMF. This replaces 
the need to use sifted white noise in complete EEMD, however, we do allow for more sophis¬ 
ticated approaches for choosing the masking signal. Note that this introduces a feedback loop 
between decomposition and demodulation. This changes EMD from simply a decomposition 
method, to a full HSA method because the lA and IF estimates of the fcth component are 
inherently computed and then used to design the masking signal for the (/c-|-l)th component. 

3. Examples using the HSA—IMF Algorithm 

In this section, we demonstrate the HSA-IMF algorithm for signal analysis and compare 
to conventional STFT using both synthetic and real-world signals. In these examples, we 
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use the proposed visualization method from Part II orthographically projected onto the 
time-frequency plane, for plotting the instantaneous parameters resulting from HSA. We 
also plot the STFT magnitude (STFTM) using the same colormap to facilitate comparisons 
[91]. The STFT is computed using a 4096-point hamming window with a 1 sample frame 
advance. We have provided matlab functions for these examples [89]. 

3.1. Synthetic Signals 

We provide two examples using synthetic signals with known underlying signal models. 
The synthetic signal examples demonstrate the proposed algorithm on a single component 
in a noiseless environment. In this case, mode mixing is not a problem and we initialize 
J = 1, a = 0.95, and /Sk = 0 in Algorithm 9. 

In the hrst example, we analyze a signal with a slow-varying AM and fast-varying FM 
given by 

(84) 

where the lA is 

a{t) = 

and the FM message is 

m{t) = 250 sin(1407rf) + 2000[exp(-4f) - 1]. (86) 

Fig. 13(a) shows the STFTM where we see classic harmonic structure resulting from the 
inherent assumption of SHCs despite (84) consisting of only a single component. Fig. 13(b) 
shows a single component in the Hilbert spectrum with a fast FM oscillation consistent 
with the underlying signal model in (84). The lA in (85), consisting of a Gaussian envelope 
centered at t = 0.5, is visible as color variation in Fig. 13(a). The FM message in (86) 
consists of a 70 Hz oscillation superimposed on a decaying exponential; this FM message 
is offset by a 3000 Hz carrier. The FM message is reflected by the vertical behavior in 
Fig. 13(b) where we see the IF sweeping from 3000 Hz to 1000 Hz with a 70 Hz oscillation. 
Due to more appropriate assumptions of an underlying component, this example shows that 
the HSA-IMF parametrization can more closely match the underlying signal model than 
traditional methods. 

In the second example, we analyze a signal with a fast-varying AM and slow-varying FM 
given by 

x(t) = jR ^a(f) exp 

where the lA is 

a{t) = A + Q sin(1007rf) -|- - sin(2007rf) (88) 

O 0 


j I lOOOvrf -|- / m{T)dT 


(87) 


x{t) = 3R < a{t) exp 


j ( OOOOvrf -|- / m{T)dT 


and the FM message is 


m(t) = 150 sin(27rf). 
45 
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Fig. 14(a) shows the STFTM where we again see classic harmonic structure resulting 
from the inherent assumption of SHCs despite (87) consisting of only a single component. 
Fig. 14(b) shows a single component in the Hilbert spectrum with a fast AM oscillation, 
captured by color variation in the plot line, consistent with the underlying signal model in 
(87). 

As described in Part II, signals composed of narrowband components can have similar 
Fourier and Hilbert spectra, however, signals composed of wideband components can have 
spectra which are quite different. Fourier analysis of wideband signals may result in multiple 
ridges in terms of spectral frequencies, which may be further decomposed into narrowband 
components. However, for many real-world signals, a wideband component, such as the AM- 
FM component in (36), is the more appropriate choice rather than multiple, narrowband 
components. The appearance of structure in the Fourier spectrum may indicate the presence 
of a wideband component (s) in the signal as demonstrated in Example 1, Fig. 13(a) (fast- 
varying IF) and in Example 2, Fig. 14(a) (fast-varying lA). 

Both decompositions, i.e. STFT and HSA-IMF are equally valid models for the real 
signal x{t )—the resulting decompositions simply correspond to the different assumptions 
of the underlying components, i.e. constant-frequency components and IMF components, 
respectively. The complex extensions assumed in STFT and implied in HSA-IMF are fun¬ 
damentally different and because of this, the quadrature signal y{t) can be different for the 
two analysis methods. This ultimately results in a different z{t) for each analysis and con¬ 
sequently different instantaneous parameterizations despite both models produce the same 
x{t) in (2). 

3.2. Real-World Signals 

We provide two examples using real-world audio signals. In Algorithm 9, we initialize 
I = 200, a = 0.95, and through experimentation, = 4 for all k. In addition, we have 
applied a 1 ms moving-average filter to smooth the IF estimate. 

In the first example, we analyze a recording (/^ = 22.050 kHz) of a single note played on 
a cello. Fig. 15(a) shows the STFTM where we again see classic harmonic structure resulting 
from the inherent assumption of SHCs. We note the fundamental frequency is approximately 
67 Hz (15 spectral lines evenly-spaced over 1000 Hz) and two dominant spectral lines at the 
second harmonic (~133 Hz) and at the hfth harmonic (~333 Hz). In addition, there is a 
brief dominant spectral line at t = 2 s corresponding to the ninth harmonic at ~600 Hz. 

Fig. 15(b) plots the five components returned by the HSA-IMF algorithm, where we see 
three dominant components. The lower two components, range in IF from ~120-140 Hz and 
from ~300-360 Hz corresponding to the dominant spectral lines in the Fourier spectrum. 
The upper component also exhibits significant energy at t = 2 at ~550-750 Hz corresponding 
to the brief dominant spectral line, i.e. ninth harmonic in the Fourier spectrum. 

In the second example, we analyze a recording (/^ = 44.1 kHz) of the word “shoot.” 
Fig. 16(a) shows the STFTM where we see the spectral energy of the fricative “SH” over 
0 < t < 0.15 s, scattered over the range 0 to 8 kHz. The spectral energy for the vowel “UW” 
over 0.15 < t < 0.25, is concentrated at a fundamental of ~230 Hz. The spectral energy for 
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the stop “T” over 0.37 < t < 0.4 s, is very weak and spread across the band and hence not 
visible in the plot. 

Fig. 16(b) plots the five components retnrned by the HSA-IMF algorithm. The “SH” 
fricative appears in three components with IF ranges of ~6000-7000 Hz (zeroth component), 
~2500-5000 Hz (first component), and ~1000-2500 Hz (second component) bnt is mostly 
captnred in the second component with qnickly varying AM and FM. The vowel “UW” is 
clearly captured in a single component near 230 Hz exhibiting some FM variation conjectured 
to be natural jitter. Unlike in the Fourier spectrum, the stop “T” is clearly captured in the 
Hilbert spectrum by the first component near t = 0.4. 

4. Discussion 

4.I. Resolving Closely-Spaced Components 

EMD has been criticized for its inability to resolve closely-spaced components and there 
have been numerous studies and analyzes on the resolving ability [93, 128, 134]. Rilling has 
investigated EMD’s (Algorithm 1) ability to resolve two tones as a function of a relative 
amplitude parameter and a relative frequency spacing parameter. This analysis describes 
regions in this parameter space where EMD returns one or two components [128]. As Rilling 
noted, the goal in EMD is not resolving closely-spaced components but rather resolving 
components that are suitably matched to an underlying signal model or compatible with 
assumptions made on the signal model [128]. 

As an example, consider two infinitely-long tones, cos (uat) and cos^uR), with Uh > tOa- 
We can express the sum of these tones as 

x(t) = 3? {exp {joJat) + exp (jubt)} (90a) 

= 3? (2 cos [(cufe - Ua) t/ 2] exp [j {ub + a;„) t/2]} . (90b) 

If Ua and Ub are sufficiently far apart, both Fourier analysis and EMD will resolve two SHCs 
as in (90a). On the other hand, if uja and ujb are not sufficiently far apart, EMD will resolve 
a single IMF as in (90b). As is well known, when uJa and oob are closely-spaced, the signal 
exhibits a beat effect. In the human auditory system, these closely-spaced tones are not 
perceived as two distinct tones but rather a single AM tone [135]. As Deering points out, 
EMD may correspond to the psychoacoustics of human hearing [93]. As Rilling points out, 
a decomposition into SHCs may not be an appropriate solution “...if the aim is to get a 
representation matched to physics (and/or perception) rather than to mathematics” [128]. 

A generalized example of this beat effect was given in Example 2 in Subsection 3.1. To 
clarify the connection to auditory perception of beating, ignore the slow-varying FM and 
hence the plot lines in Fig. 13 would be horizontal and not sinusoidal. The example, when 
analyzed with the STFT shows five closely-spaced tones as in Fig. 13(a). However, if these 
tones are closely spaced then as noted, they may be perceived as a single tone with AM 
variation. This is demonstrated with an analysis using HSA-IMF and shown in Fig. 13(b). 

A similar example regarding an FM signal, was given in Example 1 in Subsection 3.1, 
and is essentially the same waveform used in FM synthesis pioneered by drowning [62]. 
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In Chowning’s work, he expressed the FM signal as a snperposition of SHCs weighted by 
Bessel fnnctions and showed rich spectra associated with this signal. Fig. 14(a) illustrates 
this rich spectra via the presence of multiple harmonics. In the AM-FM model, such rich 
spectra may be encapsulated in a single component as illustrated in Fig. 14(b) and as has 
been suggested, a decomposition into SHCs may not be an appropriate solution to describe 
the underlying signal model. 

4-2. Remarks on Computation 

The HSA-IMF algorithm contains a triple-nested loop that incurs signihcant compu¬ 
tation depending on the signal and parameter choices. Clearly, if the signal has many 
underlying components, EMD will require more computation due to the extrema searches 
and interpolations. Unfortunately, there is no way to predict the number of components 
ahead of time. 

The outermost loop in Algorithm 9 Step 3 iteratively removes the IMF estimated by 
tone masking until termination conditions are reached. Hence, there is no way to predict 
the number of iterations required for termination of the loop. The middle loop, ensemble 
averages the IMFs returned from tone masking in Algorithm 9 Step 4. This average is 
controlled by a hxed number of trials, I (the choice of I is discussed in the next subsection). 
The inner loop results from the tone masking procedure in Algorithm 9 Step 4 calling the 
sifting algorithm twice, in Algorithm 5 Step 3. This in turn calls Algorithm 2, in which Step 
3 iteratively estimates an IMF. Within this inner loop, the step-size introduced in (79) and 
used in Step 9, also controls the speed at which termination conditions are reached. Of these 
loops, the innermost loop requires the most computation due to the search for extrema and 
interpolation. Taken together, the iterative nature of the outer and inner loops, coupled with 
computationally complex inner loop can require significant computation depending on the 
signal length and sample rate. As has been pointed out, signal oversampling is required for 
robust estimates of the IMFs further increasing computation. Finally, IMF demodulation 
in Step 5 occurs in the outermost loop and does not add significant computational burden. 
Much of this computation can occur in parallel [136, 137, 138, 139], for example, the ensemble 
averaging in Algorithm 9 Step 4 and the search for extrema in Algorithm 2 Steps 4 and 5. 

We have implemented Algorithm 9 in matlab where the trials are computed in parallel 
using a parf or loop and have timed the computation for the synthetic and real-world signal 
examples in this paper. Our PC consists of an eight-core AMD FX at 4.01 GHz with 32 
GB RAM. The results are given in Table 3 where we see for relatively short audio signals, 
decomposition may require relatively large (3 and I leading to long computation times. 

4-3. HSA-IMF Robustness 

The goal in HSA is to obtain an AM-FM decomposition with meaningful interpretation. 
Thus, proper identihcation of the underlying components is required from Algorithm 9. Two 
parameters must therefore be carefully chosen: the SNR factor flk and the number of trials 
I. As a reminder, f3k weights the additive masking signal used to mitigate mode mixing and 
I minimizes, through ensemble averaging, the influence of the additive masking signal on 
the resulting IMF; both of these parameters appear in Step 4. 
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Table 3: Benchmarks 


Example 

is (kHz) 

Duration (s) 

/3 

I 

Computation Time (s) 

Synthetic 1 

44.1 

1 

0 

1 

1.3 

Synthetic 2 

44.1 

1 

0 

1 

0.9 

Cello 

22.05 

4.21 

4 

200 

320.1 

Speech 

44.1 

0.43 

4 

200 

52.4 


In our experience with real-world signals where the underlying signal model is unknown, 
we begin by selecting = 0 and J = 1 and visualizing the resulting IMFs by plotting the 
time-real plane or time-frequency plane. Mode mixing will be evident in the time-real plane 
when the frequency of the waveform changes abruptly. Mode mixing will also be evident in 
the time-frequency plane when the IMF is similar to the illustrated IMF within the — frame 
in Fig. 11(a). If mode mixing is present, we increase I3k and I. This process is repeated 
until reasonable IMFs are obtained keeping in mind the associated computational load for 
large I. Although this process for decomposition is somewhat heuristic, such rehnement is 
typically present in all time-frequency analyzes [1], e.g. choice of window length and type 
for Fourier analysis and mother wavelet selection for wavelet analysis. 

The step-size parameter a introduced in (79) and used in Algorithm 2 Step 9, is of 
secondary importance and merely scales the trend which is removed from the signal being 
sifted. This scaling is used to minimize the impact of possible overshoot in the trend. In 
our experience, selecting a = 0.95 gives satisfactory performance, noting that lower values 
lead to additional iterations and more computation. 

The termination condition in the sifting algorithm may be too restrictively chosen hence 
preventing convergence. In our implementation we include a maximum number of iterations, 
typically 50, to guarantee an exit. As a hnal point, IMFs with excessively large values are 
omitted from the ensemble. 

5. Future Research in HSA—IMF 

In the course of this research, several avenues for further algorithm improvement are 
apparent. These include: alternative masking signals, alternate interpolators, and improve¬ 
ments to IMF demodulation. In the HSA-IMF algorithm, we use lowpass hltered noise 
as the masking signal. However, for certain signal analysis problems, more sophisticated 
masking signals have been proposed [140, 129, 130]. Ideally, the masking signals would have 
properties similar to the underlying components which would likely require prior knowledge 
of the signal model. At the heart of EMD is the iterative estimation of the IMFs which 
are completely determined by interpolation. As noted earlier, the cubic spline interpolator 
may be susceptible to overhtting which can lead to poor IMF estimates. Alternatives to cu¬ 
bic spline interpolation have been investigated including B-spline and Akima interpolators 
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[141, 142, 143, 95, 144, 145, 146], however, these do not appear to offer significant improve¬ 
ments. Nevertheless, improvements in interpolation may lead to more robust decomposition 
and demodulation. Finally, IMF demodulation requires estimation of the IF which uses 
Huang’s iterative normalization procedure. Unfortunately, the required interpolation in the 
normalization procedure can result in an overfitting of the cubic spline leading to incorrect 
instantaneous estimates. As noted earlier, changing the cubic spline interpolator in Algo¬ 
rithm 8 changes the IMF. Thus a change of the interpolator in the lA estimator requires the 
same change in the sifting algorithm. 

Although IMFs are latent AM-FM components, there are other classes of AM-FM com¬ 
ponents that are not IMFs opening up alternate possibilities. The IMF is only important 
due to its computability via sifting. Hence, it may be possible to define other useful AM-FM 
components (not defined by Cl and C2) corresponding to different assumptions on the form 
of the component (not A1 and A2), that ultimately lead to a replacement of the sifting 
algorithm as the component tracker, and thus new AM-FM decomposition methods. 

6. Summary 

In the final part of this paper, we have reported an end-to-end, solution for the estimation 
of instantaneous parameters of the AM-FM model assuming IMF components. By leverag¬ 
ing the theory developed in Part H of these papers to interpret and demodulate the results 
returned by EMD, we have provided a complete numerical method for HSA. We provided ex¬ 
amples of HSA-IMF on synthetic signals and argued that the resulting decompositions were 
more representative of the underlying signal models as compared to conventional Fourier 
analysis. Examples of HSA-IMF on real-world signals were shown to allow for alternative 
and possibly more useful interpretation of the underlying signal model. Finally, we discussed 
computational aspects of the proposed algorithm. 
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(c) 





Figure 10: This sequence of plots illustrates the steps of a first iteration of the sifting algorithm. In (a) 
the example signal composed of the components, v^o(^) ( — ) and (“); (b) the superposition of the 

components r{t) (—) and the input to the sifting algorithm; (c)-(d) the upper envelope u{t) (-•-) and lower 
envelope l{t) (-■-) of r(t); (e) average of the upper and lower envelopes e{t) « and (f) IMF estimate 

at first iteration r[t) — e{t) « 
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Figure 11: In (a) and (b), the assumed components are indicated with — and the first component or high 
frequency IMF, identified with the sifting algorithm, is indicated within the — frame. In (a), the mode 
mixing problem is apparent where we see components of disparate scales being in the same IMF (indicated 
by O) and components of similar scale in the same IMF (indicated by □). In (b), adding noise and ensemble 
averaging may assist in resolving mode mixing. 
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Figure 13: (a) STFTM and (b) Hilbert spectrum for the fast-varying FM and slow-varying AM synthetic 
signal given in (84)-(86) in Example 1 in Subsection 3.1. The wideband FM message results in harmonic 
structure under Fourier analysis and a fast-frequency-varying component under HSA. 
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Figure 14: (a) STFTM and (b) Hilbert spectrum for the fast-varying AM and slow-varying FM synthetic 
signal given in (87)-(89) in Example 2 in Subsection 3.1. The wideband lA results in harmonic structure 
under Fourier analysis and a fast-amplitude-varying component under HSA. 
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Figure 15: (a) STFTM and (b) Hilbert spectrum for the cello recording in Example 1 in Subsection 3.2. 
The lower two components in (b), range in IF from 120-140 Hz and from 300-360 Hz corresponding to the 
dominant spectral lines in the Fourier spectrum at 133 Hz and 333 Hz. The harmonics above 500 Hz in 
(a) and the upper component in (b) with IF ranging from 500-1000 Hz partially accounts for the spectral 
richness of this instrument’s note. 
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Figure 16: (a) STFTM and (b) Hilbert spectrum for the speech recording “shoot” in Example 2 in Subsection 
3.2. The “SH” fricative appears in three components with IF ranges of 6000-7000 Hz, 2500-5000 Hz, and 
1000-2500 Hz. The vowel “UW” is clearly captured in a single component near 230 Hz. Unlike in the Fourier 
spectrum, the stop “T” is clearly captured in the Hilbert spectrum by the first component near t = 0.4. 
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Conclusions 


In this paper, we began by reframing the classic complex extension problem as a Latent 
Signal Analysis (LSA) problem where the objective is to determine the complex-valued 
latent signal, z(t) from the real-valued observation, x(t) = We used this framework 

to argue against the assumption of Harmonic Correspondence (HC), and hence the use of 
Gabor’s quadrature method and the Hilbert Transform (HT) for complex extension. By 
relaxing the HC assumption, there are many choices for the quadrature and furthermore, 
many of these complex extensions can be useful in modeling physical phenomena. 

Next, we presented a theory for the Hilbert spectrum framed as generalized LSA prob¬ 
lem in which we seek a representation of z{t) consisting of a superposition of latent AM-FM 
components parameterized by a set of Instantaneous Amplitude (IA)/Instantaneous Fre¬ 
quency (IF) pairs. The use of latent AM-FM components admits many possible forms of 
the component therefore for a given x{t), there is considerable freedom in the signal model. 
We presented the analogue of the LSA problem in the frequency domain where we showed 
that a latent spectrum cannot be uniquely tied to a to the spectrum of the real observation 
because of the structure imposed by the real operator. We showed that without the HC as¬ 
sumption, there are many choices for the latent spectrum and these spectra will not have in 
Hermitian symmetry. We proposed a novel 3D visualization of the Hilbert spectrum which 
plots u(t) vs. s(t) vs. t and coloring with respect to |a(t)| and allows for the visualization of 
the instantaneous parameters. We have recast time-frequency analysis and Gabor’s analo¬ 
gies to quantum mechanics to an analysis method where uncertainty is in the quadrature 
signal and not in frequency. Further, by moving away from Simple Harmonic Components 
(SHCs) with HC, we allow for a new and powerful way to analyze nonstationary signals. 

We recognized that an Intrinsic Mode Function (IMF) is a latent AM-FM component 
and leveraged HSA theory to interpret both the IMF and Empirical Mode Decomposition 
(EMD). With this interpretation we show that the dehnition of an IMF unambiguously forces 
a unique complex extension. Furthermore, we also recognize fundamental problems with the 
Hilbert-Huang Transform (HHT), i.e. that the HT is inappropriate for IMF demodulation 
and proposed an IMF demodulation method that is compatible the the IMF dehnition. 
Finally, we utilized the EMD algorithm with our modihcations and IMF demodulation to 
calculate the lA/IF parameters of x{t), thus providing a numerical method for Hilbert 
Spectral Analysis (HSA). 
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